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ABSTRACT 


A class of impulsively started, axisymmetric, laminar jets produced by a time 
dependent point source of momentum are considered. The jets studied are different 
flows, each initially at rest in an unbounded fluid. The time dependence of the 
point source of momentum categorizes the specific type of jet under investigation. 
The study of these flows is conducted at three levels of detail discussed below. 

1. A generalized set of analytic creeping flow solutions are derived, along with 
a method of flow classification. 

2. From this generalized set, there are three specific creeping flow solutions 
which are studied in detail: the vortex ring, the round jet, and the ramp jet. This 
detailed study involves derivation of vorticity, stream function, entrainment 
diagrams, and evolution of time lines through computer animation. From the 
entrainment diagrams, critical points are derived and analyzed. It was found that 
flow geometry was dictated by the properties and location of these critical points. 
In addition, these critical points undergo bifurcation and topological transformation 
with changing Reynolds number. These bifurcations and transformations represent 
a form of transition for which specific Reynolds numbers were calculated. A state 
space trajectory was derived describing the topological behavior of these critical 
points. This state space derivation was based on continuity, and boundary 
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conditions, and was performed prior to actual solution of the momentum equation. 
From this state space derivation it was found that three states of motion are univer- 
sal for all axisymmetric jets. 

3. The third level of study examines the axisymmetric round jet which was 
solved numerically using the unsteady laminar Navier Stokes equations. These 
equations were shown to be self similar for the round jet. The boundary conditions 
used in this numerical solution are the steady solution of the round jet discovered 
by Landau (1944), and the unsteady dipole. The numerical method utilized a 
second order central difference scheme solved by an implicit matrix method. The 
matrix solver was a direct method which used a new forward-backward technique 
that greatly reduced storage requirements. The numerical method solved the round 
jet up to a Reynolds number of 30 for a 60 X 60 point mesh. From the data gen- 
erated, computer animations were produced. These animations showed each of the 
three states of motion for the round jet, including the Re = 30 case. 
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Chapter I 


INTRODUCTION 

Fascination with fluid motion has been an activity of mankind since earliest 
history. The photograph on the front page of this thesis is from the Newgrange 
passage grave in the Boyne Valley, Ireland. The tomb was constructed by preli- 
terate rule-of-thumb engineers around 2500 BC (almost as old as the Great 
Pyramid of Egypt), and is covered with swirls and spirals that might easily be 
described by the topological methods used in this thesis. One can imagine these 
ancient engineers looking into the Boyne river not far from the tomb, seeing vortices 
and eddies forming in the water, and while contemplating both their beauty and 
complexity, reproducing them on the walls of the tombs. Our work is a continua- 
tion of this ancient fascination, employing newly developed techniques for under- 
standing how and why these shapes are formed. 

1. Objectives Of The Research 

The fundamental philosophy of the research was to study flows of such simple 
geometry that they would be mathematically tractable and yet of sufficient com- 
plexity that many of the basic motions of viscous unsteady fluid flow would appear. 
This led to the study of an infinite fluid with a point momentum disturbance. 
Three simplifying assumptions that were made in studying this flow are: 
incompressibility, Newtonian fluid, and axisymmetry. 
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2. History Of The Study Of Axisymmetric Jets 

The earliest analytic realization of this type of flow was by M.J.M. Hill who 
published the spherical vortex solution of the Navier Stokes equations in 1894. The 
next major analytic breakthrough was the exact solution of steady round jet found 
by L.D. Landau in 1944 with an independent discovery by H.B. Squire in 1951. 
Probably the first numerical study of the unsteady round jet was performed by Ma 
and Ong in 1971. Their method was formulated in cylindrical coordinates and car- 
ried out with the primary aim of computing the propagation of the jet front into 
stagnant fluid. Following this work a linearized analytic solution (Stokes solution) 
of the unsteady jet was found by C. Sozou in 1979. This was preceded by a numeri- 
cal solution in spherical coordinates by C. Sozou and W.M. Pickering (1971) of the 
unsteady round jet up to a Reynolds number of 12.5. The work of Sozou and Pick- 
ering advanced understanding of the axisymmetric round jet but failed to address 
several key aspects of this flow. Most important of these was that, while the flow 
had been solved in terms of the unsteady stream function, none of the flow topology 
becomes apparent unless the flow is studied in terms of its unsteady particle paths. 
It was found by B. Cantwell in 1980 that the topology of particle paths undergoes 
critical point bifurcation and transformation at certain key Reynolds numbers. This 
discovery was made using the linearized Stokes solution found earlier by Sozou. 
The application of critical point theory to fluid mechanics was in itself not a new 
idea. Poincare himself had made applications of this theory to fluid study (circa 
1880). More recently, Oswatitsch (1958) and Lighthill (1963) classified certain criti- 
cal points which can occur near a rigid boundary. A.E. Perry and B.D. Fairlie 
(1974) were one of the earliest researchers to use critical point theory in the context 
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of A.A. Andronov’s (1971) p, q state space maps for studying fluid mechanics. 
Cantwell used the p, q map as a means of quantitatively understanding the criti- 
cal point bifurcation process as a form of transition. Questions were raised as to 
whether the technique could be applied to a numerical computation of the nonlinear 
round jet and as to the role of nonlinearity in the transition process. This is where 
the present work entered the tale of events. 

3. Problem Approach 

It was felt early in the study that a different approach from that of Sozou and 
Pickering would have to be taken. Our own goal was to fully describe the transi- 
tion process and to push the Reynolds number as high as possible based on available 
computer resources. Though our method is extremely stable it is also expensive (the 
usual tradeoff), and was stopped at Re = 30.0 although higher Reynolds numbers 
could be computed at added cost. This value is more than double the highest Rey- 
nolds number previously reported by Sozou and Pickering, which was Re = 12.50. 
Our prime objective was to apply topological methods to the numerically solved 
axisymmetric jet. As with the linearized solution, it was found that the numerical 
solution undergoes only two (topological) transitions with Reynolds number, and 
after the second and last transition (which occurs at Re = 7.54) the topology 
remains unchanged at all higher computed Reynolds numbers. It should be 
emphasized that this topological invariance is a consequence of the assumption of 
axisymmetry, and it is suspected that a non-axisymmetric round jet would be sub- 
ject to an infinite sequence of transitions and concomitant great topological com- 
plexity. As shall later be discussed, since the round jet has a similarity solution its 
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stream function or particle paths for a given Reynolds number can be represented in 
a single plot which is valid for all time. Since, at a given Reynolds number, the 
entire flow history can be very compactly stored and manipulated, it was apparent 
that this flow lent itself well to making computer animations. These computer ani- 
mations were found to be indispensable aids in understanding the flow dynamics 
and as teaching aids. 

4. Cases For Study 

Another aspect of the study was that, with the exception of autonomous flows 
like the round jet and special cases like the Hill’s spherical vortex, jet flows cannot, 
in general, be represented by self similar coordinates. On the other hand, it was 
found that all jet flows are self similar in the creeping (Re — ♦ 0) approximation 
where the Navier-Stokes equations reduce to the Stokes equations. This discovery 
essentially opened up a whole new line of inquiry in devising solutions to these 
creeping flows. In this thesis, two creeping flows were studied in considerable detail, 
along with the round jet. These flows are the vortex ring and its complementary 
flow, the ramp jet. The family of creeping flows was also studied in a wider sense 
and generalized solutions for an infinite variety a flow types are provided. The 
discovery of these generalized solutions, coupled with the capability of computer 
animation, has now provided the investigator with a new form of fluid experimen- 
tation. Though the computer can never replace the laboratory in studying fluids, 
the computer can provide very “clean” flows in the context of no outside perturba- 
tions or unwanted boundary conditions. 
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C hap ter II 

I 

I 

| PROBLEM DEFINITION 

! 

i 

! 1. Fundamental Assumptions and Equation Formulation 

The Navier Stokes equation is cast in the following form: 

I 

W = ^- + 7 +^* ( 2 . 1 ) 

where if is the velocity vector, p is pressure, p is density, f is the body 
force vector, and v is kinematic viscosity. Equation (2.1) already contains 
assumptions of incompressibility, and a Newtonian constitutive relation with 
constant viscosity. It is desired to convert Eq. (2.1) into a vorticity form so as to 
remove pressure and the body force from the equation. 

The following vector identities will be used in the development and are given 
without proof: 


!IX(VX?) = -|v (? • V) - f- V? 

(2.2) 

V x v/ = 0 

(2.3a) 

V • (v x f) = o 

(2.3b) 

v 2 ? = V (v • V) - V X(v x f) 

(2.4) 
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V X (IT X 2 ) = y(v-3-*(v-y) + (**v)?-(7-v)* (2.5) 

where jf, T are vector quantities, and / is a scalar. 

For an incompressible fluid the continuity equation becomes 

V • * = 0 . (2.6) 

Vorticity is defined as 

V X tT = o? . (2.7) 

Combining equations (2.7) and (2.3b) shows that the vorticity is solenoidal 

V • = 0 . (2.8) 

The body force vector is assumed to be derivable from a potential function 

1 = - V 7 (2-9) 

where 7 is a scalar function. Equation (2.2) is employed replacing f with H 
and combined with (2.1) via the material derivative convective term. Equation 
(2.9) is also combined with (2.1), along with (2.7) giving 

<L 1 - KX <3 = -v[ - + 7+ -~ o — - 1 +^V 2 « > - (210) 

d t Ip 2 l 

The curl is taken of (2.10). The identity of (2.3a) eliminates all of the gradient 
quantities of (2.10), and use of (2.4) with (2.6) and (2.8) proves that curl and 
the Laplacian commute for this particular problem. The result becomes 

-|y - V X (u X c3) = v v 2 5J . (2.11) 
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Equation (2.5) is now employed with (2.6) and (2.8) along with the definition of 
the material derivative, changing (2.11) into 

22- = (z3 ■ v) *+ v V 2 SJ • (2.12) 

Equation (2.12) is the basic equation for this study in vector notation. Spherical 
polar coordinates are selected as the appropriate coordinate system. Let 


u r e r + u 9 e 0 + 

(2.13) 

u r e r + u 9 e 9 + u e + 

(2.14) 


where e r , e 0 , e ^ are spherical polar unit vectors, and u/ r , u r , 

u 0 , and are scalar quantities. In spherical polar coordinates, (2.12) becomes 


dui r 

— — = ZD ^u T — It • v w r H" (2.15) 


+ v 


2 2w r 2 d {u t sm 0) 

r 2 r 2 sin 9 99 


2 dljJ p 

r 2 sin 9 9<f> 


_ 1 . 

-gj- — W • + — ( u 9 u r - ) + 


+ i' 


2 , 2 
v "» + 7 7F 




2 cos 9 


r 2 sin 2 9 r 2 sin 2 0 9<j> 


(2.16) 


9^ _ x , 

~~Qt~ = w + — [0J^u r - u^ r ) + 


(2.17) 


- 8 - 


cot 0 


{ - UfUj ) + 


+ V 


du r 




A 2 cos 0 du i 

^ j2 s i n q d<f> f2 s j n 2 q Q<t> r 2 sin 2 0 


These equations should be used in conjunction with the vector potential equation 
so as to ensure that continuity (2.6) is satisfied. The vector potential is defined 


as 


¥ = yxS . 


(2.18) 


By identity (2.3b) it can be seen that continuity is satisfied. Equation 
(2.18) is expanded into spherical polar coordinates giving 


u. = . 

' r sin 0 


1 | d{B <t> sin 6) dB e 

To d<j> 


(2.19) 


u 9 = — 


1 f 1 dB r d(rB J 


r { sin 0 d<f> 


dr 


( 2 . 20 ) 


where 


1 ( gOjfr) _ ^B r 

dr d0 


( 2 . 21 ) 


5 = B r e r + Bg e$ + B^ 


( 2 . 22 ) 


The definition of vorticity (2.7) is similarly expanded giving 
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1 

( dyu^ sin 9) 

8u e 

r sin 9 

1 89 

~d¥ 

i| 

1 du r 

d(™4,) 

r i 

sin 9 8<j> 

dr 


(2.23) 


(2.24) 


1 _ | <H™$) 8u r 

r i dr 89 


(2.25) 


Equations (2.15) to (2.25) represent the complete set of equations necessary to 
solve the unsteady nonaxisymmetric incompressible Navier Stokes problem. 
However, it is plain that these equations are very complex and a further 
simplification is desirable. This next approximation is the assumption of axisym- 
metry 


and no swirl 


8 _ 

8<f> 


= 0 , 


= 0 . 

Both of these approximations are often lumped together as the assumption of 
axisymmetry, without specifically mentioning “no swirl” which in principle could 
occur in an axisymmetric problem. With this approximation equations (2.24), 
(2.23), (2.21), (2.15), and (2.16) disappear. The remaining set becomes 


- 10 - 


~3t 


dw+ u t dwf 
Uf dr ^ r 39 


+ ^ + 


(2.26) 


cot 0 

+ u e + v 


V 2 ^ 


U! 


<t> 


r 2 sin 2 6 


r sin 9 


Ufi = — 

9 r 


3{B <j> sin 9) 
B dB 

(2.27) 

9(r BJ 

3r 

(2.28) 

3(ru $ ) 3u r | 

dr 39 1 

(2.29) 


2 . Governing Equations 

A further simplification is possible if (2.26) is combined with continuity (2.6) 
to eliminate some terms and if the Stokes stream function is defined as 


rp(r,9) = rsin9B l j ) . (2.30) 

To simplify notation, subscripts are dropped with the following redefinitions: 

W* —(jJ, U r = U, U 9 — V. 

The equations of motion become 


Ti <™> + fr (r “> + le w 


(2.31) 


i LJL 

"l r 39 


t~ — ( w sin 9 ) 

sin 9 39 


+ 
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1 dip 

i 2 sin 9 

(2.32) 

-1 dip 

r sin 9 dr 

(2.33) 

d , v du 

Tr ( ™> - m ■ 

(2.34) 


3. Self Similarity 

To bring about even further simplification it is desirable to find a self similar 
form. Self similarity can be deduced by finding a transformation for which the 
Navier Stokes equations and the boundary conditions are invariant. To find this 
transformation, (2.1) is cast into a Cartesian tensor form 


du { 

~dt 


, , ^ 
~ + ^ 


- 1 _dp_ 

p di{ 9* V dx- dx. 


(2.35) 


where x { is a Cartesian coordinate. The Navier Stokes equations can be sub- 
jected to a stretching transformation by the following change of variables: 


= 

~ f 

a U; 

(2.36) 

= 

h! 

(2.37) 

t = 

7 t' 

(2.38) 

V = 

tv' 

(2.39) 

9i = 

« 9i 

(2.40) 
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One may insert equations (2.36) - (2.40) into (2.35) giving: 


( i_) dUi> 


dt 


+ 


a 


0 


, du { 


dxj 


(2.41) 


JL 

0 


“ ~T- + )?i + 

P OX,- 


a 


13 


V 


d 2 U; 


dx\ dxj 


For (2.41) to be invariant under the transformation the following must be 


true: 


~ 2 


a a ji__ ~ a 

1 13 0 0 


(2.42) 


Replacing (2.42) into (2.36) - (2.40) forms the following Lie group: 


u, = 

x i = 
t = 

P = 
9i = 

By the Pythagorean theorem the 


t 

a a,* 

(2.43) 

~-l , 

a Xj 

(2.44) 

a t 

(2.45) 

~2 , 
a p 

(2.46) 

~ 3 , 

a 9i • 

(2.47) 

pherical 

polar radial coordinate has the 


same invariance as (2.44): 
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r = a r' . (2.48) 

Combining (2.48) with (2.43) and (2.32) the stream function transformation is 


ip = a * ip 1 . (2.49) 

The transformation of (2.43) - (2.49) shows that the Navier Stokes equation 
will admit a solution that is invariant under the stretching (2.43) - (2.47). How- 
ever not all solutions of the Navier Stokes equation have this property of invari- 
ance and to determine whether a particular solution has this property (without 
actually solving the equations) one must also examine the invariance properties of 
the boundary conditions. The boundary conditions of an unsteady axisymmetric 
jet are as follows. 

For the far field r — ► oo the solution goes to an unsteady dipole: 


_ irtR f_ s . n2 9 ( 2 . 50 ) 

4nr 

For the near field r — + 0 the solution goes to the steady Landau Squire solu- 
tion: 


ip = vr 


2 sin 2 9 
A - cos 0 


(2.51a) 


where A is a parameter which depends on the Reynolds number, Re. 


Re 2 

16?r 




A + 1 


A - 1 


(2.51b) 
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By inserting (2.48) and (2.45) into (2.50) and (2.51a), it can be seen that 
these equations satisfy the transformation of (2.49). It can now be concluded 
that an axisymmetric round jet solution can be found that is invariant under the 
stretching transformations of (2.43) - (2.47). 

The analysis just performed can also be done on the linearized version of 
(2.35) which is 


du{ 

~dt 


-1 dp , , ° *i 

"a ^ & + u H a — • 

p ox { axj dxj 


The analogue of (2.42) for the linear case is 


(2.52) 


a _ P _ ~ 

7 /? 


a 


P 


(2.53) 


The consequent group is 


«,• = a u { 


(2.54) 


= P*i 


(2.55) 




~2 

P f 


(2.56) 


a 


u I 

p = — p 

p 


(2.57) 


9i = 


a t 

TT 9i 

P 


(2.58) 
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Without the convective term the problem is underdetermined leaving two 

coefficients a and /? in the group of (2.54) - (2.58). This is a happy conse- 

quence since the additional parameter /? allows one to satisfy a wider class of 
boundary conditions. For this reason the linear (creeping) solutions are invariant 
under the two-parameter stretching transformation (2.54) - (2.58). Now that the 
invariant group is established it is possible to derive the self similar forms. Both 
(2.45), (2.44) and (2.55), (2.56) lead to: 






constant 


(2.59) 


The Pythagorean theorem with judicious selection of a constant transforms 
(2.59) into the self similar spherical polar radial coordinate 


e = 



By similar reasoning, (2.49) and (2.45) give 


(2.60) 


jL. = — — it — — = constant . (2-61) 

*2 (<rY) 7 

Again by appropriate selection of a constant, (2.61) can be converted into a 
dimensionless self similar stream function G (f, 9) 


-1 1 

= P i/ 2 f 2 G(£,0) 


^ (r, 0) 


(2.62) 


- 16 - 


where P in the case of a jet would be the strength of the momentum source. 

This process can be performed on all the dependent variables of the 
transform group. A more generalized system of self similar dimensionless vari- 
ables uses the Stokes solution group of equations, (2.54) - (2.58), where 


a 

p 2m- 1 

(2.63) 

This in turn gives the following dimensionless self similar variables 



3 1 


u (r, 9) = P 

V* f"" 2 v(t ») 

(2.64) 


3 1 

A — — _ 77 %~ — 


<2 

II 

P v 2 t 2 V(f, 9) 

(2.65) 


A -I m+ i 


ji 

II 

Pu 2 t 2 G(t,0) 

(2.66) 

u (r, 9) = 

p v- 2 r - 1 w (f, 9) . 

(2.67) 

4. The Creeping Flow, Re — ► 

0 Analysis 



The variables of (2.64) - (2.67) are equivalent to those derived from the 
Navier Stokes solution group of (2.43) - (2.47) for the case m = 0 and when 
inserted into the Navier Stokes equations eliminate the time. This shall be done 
in Chapter V where the finite difference equations will be derived. Of more 
immediate interest is the Stokes equation. The equations of (2.31) - (2.34) can be 
used to describe creeping flow or Stokes flow if the nonlinear convective terms of 
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(2,31) is removed. These equations, when combined wi^ (2.64) - (2.67), give 


8 

dZ 


2 d\v 


dt 


3 dW 


+ ; 3 — + —5— — 

* d£ sin 0 80 


sin 0 


8W 

80 


sin 2 0 


U = 


V = 


-£ W = 


£G_ 

d? 



L 

j 


m- 1)« 2 

w 

= 0 

(2.68) 

1 

8G 


(2.69) 

£ 2 sin 6 

80 


- 1 

8G 


(2.70) 

£ sin $ 

dZ 


+ a [ 1 
86 L sin 0 

8G 1 

80 J * 

(2.71) 


This system of equations can be solved, in general, and this is done in Appendix 
A. This chapter will focus only on three particular solutions which are the dipole 
form solutions for m = -1, 0, 1. The dipole form is assumed; let 


W(t,0) = sin 8J(Z) 

Equation (2.68) with (2.72) substituted becomes 


(2.72) 


— + [- + -1 — - [4 + m-ll J = 0 
d? L 2 £ J d£ J 


dt L ? 


(2.73) 


Before proceeding with the solution of (2.73) let us define the Reynolds 
number which first necessitates examining the particle path equations. 
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dr 

dt 

(2.74) 



| W 

II 

(2.75) 

Substituting (2.60), (2.64) 

and 

(2.65) into (2.74) and (2.75) gives 



d£ 

dr 

= ReU,'-| 

(2.76) 


d9_ 

dr 

= Re 2 y 

(2.77) 

where 


r = In t 

(2.78) 

and the Reynolds number is defined as 



I 



(2.79) 


With Reynolds number established we now solve (2.73) by means of the Fro- 
brenius method. The details of this procedure are lengthy and are described in 
Appendix A. The key boundary condition in all three flows is the dipole of 
(2.50). The solutions in terms of dimensionless, self similar vorticity and stream 
function along with the Reynolds number are provided for cases m = -1, 0, 1 
which correspond to the vortex ring, round jet, and ramp jet. Also provided is 
the forcing function F(t) located at the origin, which acts as a point momen- 
tum source creating the flow. F (t) is derived by recognizing that impulse is 
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i 


i 


conserved in these jet flows which provides an integral of constraint. This 
integral of constraint was examined in detail by Cantwell (1981). The precise 
mathematics for finding F(t) is found in App. A (Eqs. A-94 through A- 103). 
We can now state the three solutions for m = -1, 0, 1: 


Vortex Ring, m = -1 


F(0 = 


I) MO 


where 8 (0 = Dirac delta function, and I — impulse. 


P = — units L 4 / T 
P 


Re = 


I 

pirt 


, i 

2 


9) = ^f-7- 


-e 

sin 0 l e ~r 

16 t y/n 


G(t,e) = 


sin 2 9 

4z 


i ft £ y 1 

T Cr A n ) - -r ‘ 4 


L * 




Round Jet, m = 0 


F(t) = 


A s to 


P — — units L 4 / T 2 

P 


(2.80) 

(2.81) 

(2.82) 

(2.83) 

(2.84) 

(2.85) 

(2.86) 
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where u {t) = Heaviside step function, and J = force. 


Re = 





sin 9 

4ir 


1 e 4 


y/Z Z 


+ 


* 2 


1 - erf 





sin 2 # 

Sir 




^'^ C4 + (f‘ e ) ‘ f/ (f 


Ramp Jet, m = 1 


m 




P = — units L*/!* 

P 


Re 



2 


WU, 0) 


sin # 
2ir 






1 - erf 



(2.87) 


( 2 . 88 ) 


(2.89) 


(2.90) 


(2.91) 


(2.92) 


(2.93) 
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G (£, 9) 


sin 2 # 

4ir 



1 - erf 




£ 

2 



(2.94) 


i 

i These equations provide the basis for our analytic studies in this work. In 

passing, it should be noted that in (2.87), (the Reynolds number for the round 
jet), there is no length scale which can be derived from the governing flow 

I 

i parameters. This identifies the round jet as a very special case. It should also be 

pointed out that the round jet is exempt from the controversy over Eulerian and 
| Lagrangian frames of reference for observing flow structure. This was shown in 

i Cantwell (1981) and is a consequence of the invariance of the round jet under a 

( stretching transform as was shown earlier in this chapter. It should be pointed 

out, however, that this property does not exist for the cases of m ^ 0, except 
| in the creeping limit. It is clear now that the round jet is a very special problem 

1 which lends itself well to analysis. Its near and far field behavior are analytic 

j 

I 

and because of its self-similar nature it is exempt from many of the mathemati- 
cal difficulties encountered in most fluid studies. This makes the round jet an 

I excellent vehicle for investigating the more complicated aspects of viscous 

I 

unsteady flow. 
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Chapter III 

STATE SPACE AND CRITICAL POINTS 


In this chapter we shall develop techniques for isolating and classifying criti- 
cal points. The topology and Reynolds number behavior of these critical points 
will be examined. The critical points are derived from the particle path equa- 
tions which repeated from (2.76) through (2.79) are 


A = Re 2 U - — 
dr 2 


(3.1) 


dd_ 

dr 



Re 


A 



2 


(3.2) 


(3.3) 


1. Description Of Critical Points 

Critical points occur at coordinates 9 C such that (3.1) and (3.2) are 
both equal to zero. Equations (3.1) and (3.2) may be expanded in Taylor series 
around the critical point as 


PRECEDING PAGE BLANK NOT FftMED 
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■f- = <■(£-£«) + »(*-»«) 


(3.4) 


-f = c{(-(')+d(0-<)') 


(3.5) 


where a,b,c, and d are constants. These equations can be restructured in the 
form used by Poincare' and Andronov (1971): 


d£ 


de 


«(*-&) + b(9-$ e ) c(t-t e ) + d(e-e e ) 


= dr . 


(3.6) 


If these equations are autonomous (r does not appear explicitly except in dr), one 
may cast the equations into a second order differential equation: 


or 


where 


TF + P “T + 9 e = 
di 2 dr 


d 2 9 ^ d9 

lF + p ^ + q,> = 


0 


0 


p — — (a + d) 


(3.7) 


(3.8) 


(3.9) 


q = ad - be . 


(3.10) 


At this point the problem takes on a new perspective since it is recognized 
that (3.7) and (3.8) are both equations for the damped harmonic oscillator where 




jf fi * ” : . 
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the damping and spring force would determine p and q (which is the same for 
both equations). The parameters p and q can be found by using (3.9) and 
(3.10) and derivatives of (3.1) and (3.2), i.e., 


a 


b 


c 


d 


d 


d_ 

d$ 


d_ 

dt 


d 

dd 


A 

dr 




A 

dr 






(3.11) 


(3.12) 


(3.13) 


(3.14) 


Equations (3.7) and (3.8) have the same characteristic equation 


~ 2 

a + p a + q = 0. 


The solution of this quadratic equation is 


(3.15) 


a 


= - J ( P±\V- 4?) 


(3.16) 


Using (3.16) we can divide up p , q space into particular regions which will 



- 26- 


have distinctive topologies. Figure 3-1 shows these particular regions in p , q 
space with their critical points displayed. This plot shown in Fig. 3-1 goes by 
many different names, i.e., p , q space, stability diagram, Poincare' map, etc. 

The four critical points of prime interest in this paper correspond to values 
of p > 0. The first of these is the stable focus which exists for q > p 2 / 4. The 
stable focus represents a region in p , q space where (3.16) is imaginary. 

Next is the star point for p 2 = 4 q where the radical of (3.16) goes to zero. 
The stable node is for p 2 /4 > q > 0, which represent a real value of (3.16). 
Finally, the saddle point is for q < 0. 

As an example of how one of these topologies is formed we will examine the 
case of the stable focus. In this case the solutions of the two equations (3.7) and 
(3.8) are behaving as decaying sinusoidal functions, both equations with the 
same frequency and rate of decay. These equations can be coupled together to 
form a Lissajou figure where two coordinates decaying as damped sinusoids form 
a two dimensional spiral, spiraling inward with time. In the case of the stable 
center (p = 0) solution (3.7) and (3.8) are not decaying, but are pure sinusoids. 
In this case one has a family of concentric closed curves (circles or ellipses). 

2. Universality Of The p, q Trajectory 

The location and identification of critical points of solution is only part of 
the analysis. As shall soon be demonstrated, it is possible to predict the flow 


UNSTABLE FOCUS 
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i 







Figure 3-1 Critical Point Topologies in p, q Space 
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trajectory in p, q space prior to solution of the momentum equation. Having 
such a p, q trajectory, one can anticipate the flow topology and the manner in 
which it will change with Reynolds number without having to solve difficult 
equations. This method is developed by first combining (3.1) and (3.2) with (3.4) 
and (3.5) 


U 


1 

Re 2 


«<f - ? e > + -»«) + { 


(3.17) 


V = [«({-&) + (3.18) 

where a, b, c, and d are as yet undetermined. 

Continuity is the principal equation needed for this analysis. In spherical 
polar coordinates, continuity is 


I 

r 


d_ 

dr 


(ru) + 


1 

sin 8 


m {v sin e) 


0 . 


(3.19) 


Equation (3.19) is easily converted to nondimensional, self similar form by using 
the earlier velocity definitions, giving 


2V + (™ + 

d( 


dV . cos 0 


de 


+ 


sin 8 


= 0 


V 


(3.20) 
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We can now take (3.17) and (3.18) with appropriate derivatives and combine 
them with (3.20) giving 

sin 0 1 2s (( - (,) + 26 (0 - O + ?(| + a + <j) J (3.21) 
+ { cos «[<«-{,) + .<*-*«)] = 0 . 

In examining the p, q behavior of critical points at arbitrary locations in the 
physical space of the flow it is sufficient to study three possible cases in the upper 
half plane (which is effectively the whole plane, since the problem is symmetric 
by definition), namely 0 = 0, 0 < 0 < x, and 0 = ir for any £. 

Case 1: 6 C = 0 

cos 6 ~ 1 — — 

2 

sin 0 9 . 

The boundary condition due to axisymmetry is: V = 0 therefore by (3.18) 
c = 0. Inserting the above into (3.21) 

2a(^-^) + 260 + e(| + a+</) = 0 • 
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Now let 6 — * 9 C , £ — ► £ e , utilize (3.9), the above equation becomes: 

P = d + “ • (3.22) 

Now combining equations (3.9), (3.10), (3.22), and keeping in mind that 

c = 0 for this case, we find that 


<7 = 



(3.23) 


Case 2: 0 < 0 C < tt 

For this case sin 0 is always nonzero. Let 6 — ► 0 e , f in (3.21) 

which becomes: 


P = | • (3.24) 

Case 3: 6 C = ir 

cos 9 tt -1 + ~ ^ 

2 

sin 9 tt ir - 0 


Boundary condition due to axisymmetry: V = 0, therefore by (3.18) c = 0. 
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Insert the above in (3.21) and one obtains: 

2a (f - Zc) + 26 (0 “ *) + £ [ y + a + rf) + f | 1 - ^ ~ ^ j = 0 • 

Now let 8 —*■ 9 e , f — < ► £ e and utilize (3.11) and the above equation becomes 

V = d + | (3.25) 

which is the identical result found for the case of 9 e = 0. Likewise the same 
equation as (3.23) will be generated for 0 e = ir. This should be no major 
surprise since nowhere in this p, q theory has a direction-of-flow been assumed, 
so 9 e = 0, 9 C = 7 T should have identical trajectories. The trajectories 
derived from (3.23) and (3.24) can now be plotted, see Figure 3-2. 

We are now in a position to make some rather sweeping statements about 
the flow. For example: The class of axisymmetric flows under study cannot form 
unstable nodes, unstable foci, or stable centers. A stable focus will only form 
off the axis of the flow (p = 3/2). Likewise, a star point can form off-axis at 
( p , q ) = (3/2, 9/16) and on axis at (1, 1/4). The latter point corresponds 
to zero flow. These statements can be made before the momentum equation is 
considered. The results come solely from continuity and similarity which is a 
consequence of the nature of the time dependent force which drives the flow. 



Figure 3-2. State Space With Axisymmetric Trajectory 
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Figure 3-2 is valid for all self-similar axisymmetric flows. 

However, for a specific flow such as the round jet there is some extraneous 
information that will have to be removed. (Just as there is in a generalized solu- 
tion to a differential equation prior to application of boundary conditions.) In 
determining the relevant parts of Figure 3-2 it is recognized that the geometry of 
a round jet is bracketed by two limiting aspects: The near field aspect, which is 
the steady Landau-Squire solution (2.51a), and the far field aspect, which is the 
unsteady dipole (2.50). Each of these solutions can be cast into the particle path 
form of (3.1) and (3.2) through (2.69) and (2.70), these equations in turn can be 
operated on by (3.11) - (3.14) which provide us p and q through (3.9) and 
(3.10). The results are the following. 


Landau-Squire Solution: 


where A is defined in (2.51b), 


1 2 £ sin 2 6 

Re 2 A — cos 6 

p = 5/4, and q = 1/4. 


(3.26) 


Unsteady Dipole: 


Coo(^) = 


sin 2 0 
4tt£ 


(3.27) 


where p = 7/4, and q = -1/2. See also Cantwell (1981, pp. 378-379). Both of 
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these points are on the parabola given by (3.23) which describes critical points 
on the axis of the flow. 

3. Topological Transition 

In getting from (p , q ) = (5/4, 1/4 ) to (7/4, -1/2) the trajectory 
must pass through (3/2, 0) where the line p = 3/2 for the case of “off-axis 
critical points” intersects the parabola for “on-axis critical points.” Now the 
question arises whether the line p = 3/2 can have negative or positive values 
of q for the round jet. 

It has been found in other studies on critical points that the total number of 
node points subtracted from the total number of saddle points is an invariant of 
the flow. If the p , q trajectory starts out at (5/4, 1/4) with a stable node 
and ends at (5/4, 1/4) with a saddle point, then somewhere in route two node 
points will have to be produced so that the total remains unchanged. The inter- 
section at (3/2, 0) provides such an opportunity since just as the on-axis para- 
bola crosses from the node domain to the saddle domain in Fig. 3-2, it can pro- 
duce through bifurcation new stable nodes on the p =3/2 line. However, 
one should recognize that the two off-axis stable nodes actually represent a single 
critical line centered on the axis of the jet. When viewed from a perspective 
looking along the axis of the jet, the off-axis critical line is a ring whose parame- 
ter increases with Reynolds number. So the process of bifurcation involves the 
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splitting of a single on-axis stable node critical point into an off-axis stable nodal 
line. This change in the topology of the flow represents a form of transition since 
the changes occur at specific values of the Reynolds number. 

If we redraw Fig. 3-2 for the specific case of the round jet, see Fig. 3-3, we 
see that at q = 9/16 the off-axis critical point (line) becomes a star point and 
then for q > 9/16 the off-axis critical point becomes a stable focus. This topo- 
logical transformation at (p, q) = (3/2, 9/16) represents a second transition 
after the first transition already mentioned. 

4. The Three Flow States For The Round Jet 

It is now recognized that the round jet has three possible topologies or states 
which are partitioned by two transition Reynolds numbers Rej and Re 2 . Let 
us describe these states: 

State 1: 0 < Re < Re x 

One on-axis stable node. 

State 2: Re j < Re < Re 2 

One on-axis saddle point and, one off-axis stable node crit- 
ical line forming a ring around the axis of symmetry. 


State 3: 


Re 2 ^ Re 
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One on-axis saddle point and, one off-axis stable focus crit- 
ical line forming a ring around the axis of symmetry. 

As the reader can now see, the method of p, q analysis is very powerful, for 
without having defined or solved a momentum equation we know all the possible 
topologies of viscous axisymmetric round jet and all the possible modes of topo- 
logical transition. 

Before closing this chapter we should examine a paradox that arises from 
this method. The state of rest (zero flow) is described by the on-axis star at (p, 
q) = (1, 1/4). This may not be immediately apparent to the reader unless it is 
recalled that the coordinate system from which the p, q trace is derived ( f, 9 ) 
has time embedded within it (see Eq. 3.1). Therefore if one were to draw a circle 
in physical space this same circle would appear to be shrinking in (f, 0 ) space 
(like drawing a circle on an inflated balloon and letting the air leak out). 

It so happens that only a star point has this feature of the flow field rushing 
in a purely radial direction towards the critical point in the center. Unfor- 
tunately the p , q trace displayed in Fig. 3-3 never reaches the star point at 
(p, q) = (1, 1/4). The only star point admitted by an axisymmetric jet is 
the star critical line at (3/2, 9/16) which is off-axis and not representative of 
a stationary flow. The closest an axisymmetric jet gets to (1, 1/4) is the 
Landau-Squire solution of (5/4, 1/4). What compounds the paradox is that 
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Figure 3-3. State Space With Round Jet Trajectory. 
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the Landau-Squire solution admits the zero flow case, Re = 0, A — ► oo but its 
p , q is still (5/4, 1/4). In other -words the limit Re — ► 0 corresponds to 
a topologically different flow from the case Re = 0. 

Lastly, we should address the question of nonautonomous flows, that is, 
flows where the Reynolds number is time dependent. The p , q method of 
analysis can be used in a different manner on nonautonomous flows. There is 
nothing that prevents one from freezing the particle paths for a given time, utiliz- 
ing equations (3.11) - (3.14), and calculating an instantaneous p , q. If one 
presupposes, as an assumption, that the self similar coordinate is valid for nonau- 
tonomous flows of all Reynolds number (which is true in general only as Re — ► 
0), then a trajectory identical to Fig. 3-2 can be produced. 

The difficulty with nonautonomous problems is that the particle path equa- 
tions are constantly changing. In both the autonomous and nonautonomous 
problem one can take the particle path equations (3.1) and (3.2) and generate a 
field of vectors of unit length which would show the direction of particle paths at 
a given point and at a given instant, tunnel model). A plot of these vectors is 
called an “entrainment diagram”. The entrainment diagram in the autonomous 
problem does not change with time. Therefore in the autonomous case if one 
takes a particular point and follows its trajectory in time, its trajectory will 
always be tangent to the vectors of its entrainment diagram. 
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This trajectory of a particular point through time will be referred to as an 
“integrated particle path”. However, in the nonautonomous problem the entrain- 
ment diagram is constantly changing. One can take the entrainment diagram of 
time and overlay it with a particle path of 0 < t x < t 2 where times 0 
and t 2 mark the end points of the particle path trajectory. The vectors on the 
i entrainment diagram need not be tangent to the integrated particle path trajec- 

i tory except at the point of the particle path trajectory that was made at 

t v Each of these instantaneous entrainment diagrams will have a flow topology 
, with critical points that can be studied in a p , q context. 

| However, this p , q context is different from the autonomous example 

I 

| because time can play a part. One must recast (3.4) and (3.5) as 

= a (t - (') + W - O') (3.28) 

t d a 

= '((-(') + * (0-<y (3.29) 

da 

where a is just a dummy independent variable. One can produce a “pseudo 
particle path” for an instantaneous time by integrating (3.28) and (3.29) in a . 
The “pseudo particle paths” have nothing to do with real particle paths but 
merely serve as a vehicle for observing instantaneous topologies. Standard p , q 
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analysis can be applied to this “pseudo particle path” and it is from here that 
instantaneous nonautonomous p , q values have their basis. 
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Chapter IV 

LINEARIZED RESULTS 

1. Creeping Flow Particle Path Equations 

In Chapter II the linearized solutions in terms of vorticity and stream func- 
tion are developed for the cases of the vortex ring, round jet, and ramp jet. In 
Chapter III the chief methods of analysis were developed in terms of the particle 
path and p , q plot. As the reader is probably aware, in unsteady fluid 
mechanics particle paths are the most desirable means for displaying fluid 
motion. All three of the linear flows appear as simple dipoles in the stream func- 
tion and vorticity. In terms of these variables much of the flow topology and 
influence of critical points is not apparent. We shall use Eqs. (2.69), (2.70), 
(2.76) and (2.77) so that the particle path equations become 


d§_ = 

Re 2 

3G e 

(4.1) 

dr 

f 2 sin 9 

39 2 

d0_ _ 

Re 2 

_ 3G_ 

(4.2) 

dr 

f 2 sin 

9 


where r = In t. The stream functions derived in Chapter II Eqs. (2.84), 
(2.89), and (2.94) are inserted into (4.1) and (4.2) giving: 


Vortex Ring (m = -1): F 


d£ _ Re 2 cos 9 
dr 2x£ 2 


dd _ Re 2 sin 9 

dr 4 ;r£ 3 


Round Jet (m = 0): R 


d£ Re 2 cos 0 

dr ~ 2x£ 2 


dd Re 2 sin 9 

dr ~ 4 ir? 


Ramp Jet (m = 1): F 
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It should be reemphasized that the Reynolds number functional form with 
respect to time is different for each flow (see Eq. 3.3). 

2. Critical Point Analysis Of The Creeping Solutions 

Equations (4.3) through (4.8) represent the actual equations used to plot par- 
ticle paths and locate and characterize critical points. These are integrated 
numerically to produce computer animations of the respective flows. Critical 
points are found by setting the right-hand sides of the particle path (4.1) and 
(4.2) equal to zero, and finding the value of (f e , 9 e ) for a given Reynolds 
number. The method of analysis is the same for all three flows. The vortex ring 
will be used as an example case since it is the most mathematically compact. If 
(4.4) is examined, one finds that it will go to zero for 0 = 0 or where the 
expression in parenthesis of (4.4) goes to zero. This expression can be cast as a 
transcendental equation for £ c . 
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= 



(4.9) 


This equation is in an iteratively stable form and quickly converges to a root 

A. 

£ e = 3.0224, (which shall be called £ e ). Where (4.4) goes to zero, so the root 
which causes (4.3) to go to zero must also be found to isolate the critical point. 
Setting the left-hand side of (4.3) to zero gives 


Re 2 





i 

v/?r 


dl 

, 4 


cos 0, 


(4.10) 


For the on- axis case of 0 C = 0, (4.10) becomes a straightforward function 

A 

R e = /(£ c )> unt il £e == £ e • At this point bifurcation is possible since (4.4) is 
zero for 0 C ^ 0. 

The bifurcation Reynolds number can be calculated by setting £ e in (4.10) 

A 

equal to f e . Thus 


0 e = ± arc cos 


18.1749 | 2 
Re ) 


(4.11) 


It was found that for 6 e = 0, (which is the angle at bifurcation) the Reynolds 
number is Re = 18.1749. This is the first transition Reynolds number (Rej as 
referred to in Chapter HI). Equation (4.11) also establishes the precise location of 
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the off-axis critical point after bifurcation since it provides 9 e for any given 

A 

Re and it is known that £ e — £ ,. for the off-axis critical point. 

For the on-axis critical point for Re > Rej, we continue to use (4.10) 

A 

with 0 e = 0, except now £ e >£ e . The method of determining the location 
of critical points is the same for the other two flows as just shown for the vortex 
ring. Equation (4.11) can be generalized as 


0 e 


± arc cos 


Re ^ 2 
Re ) 


(4.12) 


The results of this method when performed on all three flows is shown in Table 
4-1. 


Table 4-1 

First Transition Constants 


Flow Type 

M 

A. 

Rej 

Vortex Ring 

-1 

3.0224 

18.1749 

Round Jet 

0 

1.7633 

6.7806 

Ramp Jet 

1 

1.2821 

3.7386 
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One difficulty that is immediately apparent in Table 4-1 is that the values of 
Rej are beyond the domain of validity for a creeping flow approximation. In 
the case of the round jet this is not as bad as it might seem since in Chapter II it 
was shown that the topology for autonomous flows is a consequence of continuity. 
As shall later be shown, the main consequence of the creeping approximation in 
the round jet geometry is that the jet is of shorter length along the axis of sym- 
metry, the spreading angle of the off-axis critical points is too large, and the 
value of Rej, (and Re 2 ) is different. However, the basic topology in this 
approximation is correct. In the case of the nonautonomous flows (ramp jet and 
vortex ring), there are more significant problems. As mentioned in Chapter II 
these flows are self similar in f only in the creeping approximation and as the 
creeping approximation loses validity so does the use of this coordinate and the 
topologies that are a consequence. This fact has been observed by Glezer (1982) 
in his experimental work for the turbulent vortex ring. They used the self simi- 
lar variable of r/[It j 1 / 4 , (where / is the impulse), and also found an additional 
critical point on the axis of symmetry. Thus the nonautonomous creeping solu- 
tion provides only a partially complete topology and transition behavior. 

One other aspect of Table 4-1 is that Rej seems to be diminishing as m 
increases. Though it has not yet been determined, it might be found that for 
sufficiently large m, Re! will occur for values that are within the creeping 
approximation which would represent very interesting subjects for low Reynolds 
number experiments. 
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We now return to our example problem of the vortex ring and examine the 
situation for Re > Re^ The values of p , q can be calculated for the vortex 
ring by taking (4.3) and (4.4), inserting them into (3.13) - (3.16) and utilizing 
the definitions of (3.11) and (3.12). The resultant equations require considerable 
algebraic manipulation and are quite long. However, they can be simplified by 
recognizing that in the on-axis case, {9 e = 0), q is given by (3.23) which is 

solely a function of p. Therefore only p needs to be calculated. This is 

given by 


Ptf c =0 


_3 _ _Re^_ 
2 " 4xg 



(4.13) 


As was shown earlier, is solely a function of Re for 6 e = 0 and can be 

found by using (4.10). Next, for 0 e ^ 0 it is already known that p = 3/2, 

/ 

so only q needs to be found. The very complicated equation for q is reduced 


by recalling that = £ e for 6 e ^ 0, and utilizing (4.12). After much alge- 
bra and lengthy hand computations one finds that 


Qa = F + A Re 4 

where for the vortex ring: 


(4.14) 


r = -0.3209 5 

= 2.9413 X 10" 6 


A 
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Equation (4.14) is the basis for calculating the second transition where the off-axis 
critical point transforms from a stable node to a stable focus. This transition 
occurs at the parabola shown in Fig. 3-2, which is 

1 = -f- • (4.15) 

One can now set p = 3/2 yielding q = 9/16 and insert this result into 

(4.14) giving Reo = 23.4105, which is the second transition Reynolds number 
for the vortex ring. The angle 9 e at which the second transition occurs can be 
calculated from (4.12). The method used for the vortex ring can be applied to 
the round jet and ramp jet -- the results of which are shown in Table 4-2, (Rej 
is restated from Table 4-1 for comparison). T and A in Table 4-2 refer to Eq. 

(4.14) . 


Table 4-2 

Second Transition Constants 


Flow Type 

M 

Rej 

Re 2 

r 

A 

0 C @ Rea 
(deg) 

Vortex Ring 

-1 

18.1749 

23.4105 

-0.32095 

2.9413 X 10" 6 

52.9343 

Round Jet 

0 

6.7804 

10.0909 

-0.14405 

6.8143 X 10 5 

63.1605 

Ramp Jet 

1 

3.7386 

5.7887 

-0.11849 

6.0652 X 10" 4 

65.3468 
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For Re > Re 2 there are no further topological changes for the axisym- 

metric problem. As Re -* oo the off-axis critical point remains a stable focus 

and its q value goes to infinity as can be seen in (4.14). The off-axis critical 
point angle 9 C goes to n/2 for infinite Reynolds number as is shown by 
(4.12). This, of course, is a consequence of the creeping approximation. For the 
limiting behavior of the on-axis critical point p , q values, one uses (4.13) with 
Re eliminated by substituting (4.10). The resultant equation will yield values 
for p in the limit of Re — ► 0 where —*■ 0, and Re — ► oo where 
£ e -* oo. This same procedure can be used in the round jet and ramp jet. 
The q values can be determined in all three flows by employing (3.23) and 
inserting the respective values of p. The result of this analysis is that for the 
vortex ring limit Re — *• 0, (p , q ) = (1, 1/4), and in both the round jet 
and ramp jet limit Re — ► 0, (p , q ) = (5/4, 1/4) is found. For the case 

of Re — ► oo in all three jets (p , q ) = (7/4, -1/2). This is consistent 

with intuition which would suggest a far-field dipole behavior for all three flows. 
Now p , q plots can be drawn for all three flows. These are shown in Figure 
4-1. With the p , q plots of Fig. 4-1 and the equations developed in this 
chapter, one can describe how the respective flows will behave without integrat- 
ing the particle path equations. The vortex ring will have an infinite Reynolds 
number at t = 0 (recall that Re is proportional to 1/f 1 / 2 ). Its off-axis critical 
point will have an infinite q value so the flow will be rolling up very rapidly 


around the stable focus. 
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Fig. 4-1: State Space Showing p (Re), q (Re) Trajectories of a) Vortex 

Ring, b) Round Jet, c) Ramp Jet. 
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The critical point will be at 9 e = 90° initially but its angle will diminish 
as time progresses just as q decreases. As the q value approaches the star 
point on the q = p 2/4 parabola, the rate of roll-up will be so slow as to be 
imperceptible. When (p , q ) goes below (3/2, 1/4) and the second transition 
occurs (actually the first transition chronologically in the vortex ring), roll-up 
will cease altogether, only straining and translation occurs thereafter. As time 
progresses, 0 e for the off-axis point diminishes until Re = Re t is reached. 
After this “first” transition, the now simplified topology with one on-axis critical 
point revolves with increasing q but diminishing p along the p , q parabola 
for the on-axis critical point. Eventually after infinite passage of time the (p, q) 
value of (1, 1/4) is attained (the zero flow star) and the flow comes to rest. In 
the round jet the p , q trajectory and critical points are static because the 
flow is autonomous. This permits the realization of flows which stay in the same 
topological state for all time. Thus it is possible to study a state of motion which 
is only a transient phenomenon in a nonautonomous flow. The ramp jet is the 
complement of the vortex ring (due to Re being proportional to f 1 / 2 ), as com- 
pared to r 1 / 2 except for the fact that like the round jet this flow also has the 
start paradox at (p , q ) = (5/4, 1/4). The single stable node later bifurcates at 
the first transition (which is now chronologically correct). The on-axis saddle 
goes to (7/4, -1/2) at infinite time. The off- axis stable node point goes 
through second transition to a stable focus and attains infinite q at 6 e — 90 ° 
at infinite time. With the ramp jet we can anticipate an initially sluggish flow 



- 52 - 


starting to roll up with the rate of roll-up getting faster with time (the opposite 
of the vortex ring). 
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Chapter V 


THE NUMERICAL METHOD 


1. Derivation Of The Finite Difference Equation. 

In the previous chapters and in Appendix A the Stokes and Navier Stokes 
equations were examined in the context of analytic solutions. It was shown in 
Chapter II that the round jet (m = 0) has a self similar solution to the Navier 
Stokes equation. The Navier Stokes equation and vorticity definition can be 
recast into this self similar form by inserting equations (2.66) and (2.67) into 
(2.31) - (2.34) using the self similar coordinate of (2.60) giving 


C 2 


d 2 W , 

ae 



, i! 
86 2 


aw 

dt 


+ 


c2 _ i , i !£ 

^ sin 2 6 Z sin 6 dd 


cos e del d 2 W 

■ \\ -i_ 

sin 2 6 &Z 80 2 


+ 


cot 6 + 


1 

sin 9 


dG 

dZ 


dW 
d 6 


(5.1) 


W = 


-1 

f 3 sin 6 



+ 


fG_ 
d 9 2 


+ cot 9 


dG 

d9 


(5.2) 


Experience has shown that (5.1) and (5.2) need to be cast into a more gen- 



eralized form before conversion to a finite difference equation. The more general- 


ized form for (5.1) is 


* Hf*«' 

- MS* [■«]#- • « 

where 

M = sin 2 9 - 2£ 2 sin 2 6 - £/? sin 9 -^7 ~r (5.4 

2 da do 

eZ 

N = ±-a sin 2 0 (5.5 


B — — 2 f 3 sin 2 9 + a M 

2h ( o£ 


(5.6 


C 


d 2 a . da 


t 6 sin“0 —— + M + a 


df 2 d£ 


2£ sin 2 0 + sin 9 ^ f 
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d 2 S . 85 | , . „ „ . . A 80 „ . a 8Z ) 1 \ 

= if 1 sur '' w M [ 1 l * ra ' cos ' T sm ’ Ilf " T p ' w ) J / 


The dependent variables used are 


(5.9) 


m,») = 4 »(« w) 

f 2 

G(f,») - /MOW) 


(5.10) 


(5.11) 



— i 


II 

a 

+i - er/ i 

— "i 

Uy>| <M 


(5.12) 




zl 1 



II 

£ 

2 

•Jx + 

[I.il 
( 2 1 

er/ U) 


(5.13) 


The equation (5.12) is the radial component of the creeping vorticity solution of 

(2.88) . Equation (5.13) is the radial component of the creeping stream function of 

(2.89) . Both of these relations will be referred to as “nonlinear scale functions”. 
The singularity of 1/f in (5.13) is removable as f —*■ 0 because for small 

e 


erf 


1 

2 


v/tt 


The coordinate 6 is an angular coordinate for use in mesh stretching 


based on 0 and the Reynolds number 
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6 = 6 (9, Be) . 


(5.14) 


Equation (5.14) was introduced into (5.1) by the following equations 


OW 

09 


08_ 0W_ 
09 08 


(5.15) 


cPW 

oft 


08 

09 


0 2 W jft6_ OW 

08 2 Oft 98 


(5.16) 


The specific “theta stretching algorithm” will be discussed later in this chapter. 
The more generalized form of the vorticity definition (5.2) is 


[j?vv] 0 + [2*5] ff + [5] z+[i?d] 
+ § = * 


(5.17) 


N = I P sin 9 


(5.18) 


~ 1 . a -0 Op 

B = T sm«f — 


(5.19) 


C = ft sin 9 


Or p 

Oft 


(5.20) 


D = 


P . A 06 
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(5.21) 


E = 


2k 
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sin 9 — — - cos 9 

Oft 


06 

OO 


(5.22) 
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R = -^aTsin 2 # . (5.23) 

The variables h and k which appear in (5.3) - (5.22) are at this stage just 
dummy variables. If we subject (5.3) and (5.17) to a second-order finite difference 
we have the following 


(N+B)Y i+lJ + (N-B)Y i _ lJ +[C-2(N+D)]Y i<J . 

+ (£ + E)Y ij+l + (D- E)Y { j_ t — 0 ; (5.25) 


1 


" 


■ 

A# 

N+B 

Z i+ U + 

A# 

N-B 

+ 

A< A« «W 

C - 2(N + D ) 


D + E 


hj+x + l*> ~ E 1 = R 


(5.26) 


Equations (5.25) and (5.26) form the basis of a “two step” method of solution 
where vorticity is solved for first and then used to solve for the stream function. 
Terms £ and 8 will now be defined in a finite difference context 


Z — ih (5.27) 

6 = jk (5.28) 

where » and j are integers marking a particular node in an n X n square mesh. 
We can see now that h, and k are step sizes in the f and 6 direction. The 

coefficients N, B, C, etc., are the same as defined for the partial differential 

equations in (5.4) - (5.9). The following derivatives have a finite difference 
approximation 
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dZ 

z i+ Ij “ Z i-hi 

(5.29) 


2h 

dZ 

Z i,j+l - Z iJ - 1 

(5.30) 

08 

~ 2k 


Equation (5.29) is employed in (5.7) and (5.9) while (5.30) is used in (5.4) and 
(5.7). The equations (5.25) and (5.26) represent the actual finite difference equa- 
tions used in the software. 

2. Boundary Conditions In The Near Field 

At this point we should examine the boundary conditions in the context of 
this numerical scheme. Figure 5-1 shows the computational mesh which is n 
points by n points, with the boundary conditions as shown in Figure 5-1. The 
computational domain is actually a semicircle. 

However, if one were to visualize this semicircle as a spread-out Japanese fan 
one could imagine taking the pivot pin out of the base of the fan and unfolding it 
into a square piece of paper. This is represented in Fig. 5-1 where the radial and 
angular component have been mapped into a Cartesian system with the line 
£ = o corresponding to the point momentum source. The boundaries at 
0 = 0 and 0 = x are both on the axis of the jet and have the same boundary 
conditions of vorticity and stream function being equal to zero. For the near 
field ( £ = 0) , we employ the Landau-Squire solution. 



NEAR FIELD LANDAU-SQUIRE BOUNDARY 



FAR FIELD MULTIPOLE BOUNDARY 
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The dimensionless self similar stream function of the Landau-Squire solution 
is provided by (3.26) which when combined with (5.13) and (5.11) gives 


limit Z(£,6) 
*-o 


16 7T sin 2 0 
Re 2 ( A - cos 0) 


(5.31) 


where A is defined by (2.51b). 

For the vorticity at £ — < ► 0, we start with (2.51a) and insert this into 
(2.34) via (2.32) - (2.33) and then develop an expression for vorticity in physi- 
cal coordinates (Eq. A-115 of App. A). This result is cast into a dimensionless 
self similar form via (2.67). When combined with (5.10) and (5.12), one obtains 


limit T (£,<5) 
* - o 


16 7T (A 2 - 1) sin 0 
Re 2 (A - cos 0) 3 


(5.32) 


Equations (5.31) and (5.32) are the f -*• 0 boundary conditions used in the 
software although (5.32) needs to be algebraically manipulated to show its 
equivalence. It is convenient to tie up loose ends at this point regarding the 
parameter A of the Landau-Squire solution which is related to the Reynolds 
number in (2.51b). In Fig. 5-2 we see that as Re -*• 0, A -*• oo and as 
Re -* oo, A -* 1. Equation (2.51b) can be expanded into an infinite series 
giving, 


Re 2 
16 JT 


E 


71=0 


9 + 8n 1 
9 + 6n A 2n+l 


(5.33) 


F rom (5.33) we can see the asymptotic expansion in (5.34) 
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Figure 5-2. Plot of A vs Re. 
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Re 2 
16 n 


lim -7- . 

A -00 A 


(534) 


If (5.34) is inserted in (5.31) or (5.32) and A — ► 00 it can be observed 
that the Reynolds number dependence cancels out in both equations as Re — ► 0. 
While it is straightforward to calculate Re for a given A it is remarkably 
difficult to calculate A for given Re. No stable implicit iterating scheme has 
been discovered for (2.51b). Because (2.51b) is basically a corner function, one 
finds that Newton-Raphson schemes and bisection schemes will fail unless one is 
very careful. This is due to the fact that for one leg of the function a slight 
change in A causes an extreme change in Re while for the other leg a slight 

change in Re causes an extreme change in A (see Fig. 5-2). Therefore if an 

interating scheme is tuned for one leg it will fail for the other and no scheme 
seems to work well for both. The approach eventually used was to employ a 
brute force bisection iterator for small A and to employ the asymptotic relation 
(5.34) for large A (A > 50.2880, Re < 1.0). Though this method works it is 

slow and inefficient. The author feels that a better method can be found and 

suggests either two iterators customized for Re > 4.5 and Re < 4.5 (which is 
roughly the turning point in this corner function) or attacking (5.33) by using 
the relation A = 1 + e where e is assumed to be small (efforts at this 
approach have so far been fruitless). 


3. Boundary Condition In The Far Field 

The remaining boundary condition that needs to be examined is the far field 
(r — * co) boundary. Placing the dipole boundary as expressed by (2.50) at a 
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finite location was tried initially but undesirable oscillations and degraded conver- 
gence rates were a consequence. One positive aspect that was noted though, was 
that the near-field behavior of the solution including critical point location q 
values, etc. was strongly dominated by the Landau-Squire solution boundary of 
(5.32). One could have almost anything as a far field boundary (including gar- 
bage as was discovered) and provided it was far enough away, have little effect 
on the near field solution. This aspect emboldened us to try different approaches 
to the far field boundary. The far field vorticity was quite straightforward and 
was simply set to zero. Vorticity dies off as a variation of some Gaussian func- 
tion, in an unbounded problem, and has been observed to drop by tens of orders 
of magnitude in the range of f used. 

The approach we ultimately decided to use for the far field boundary condi- 
tion involved the axisymmetric irrotational multipole solution described by Eq. 
(A- 113). The far field portion of (A-113) was selected, combined with (2.66) to 
provide the self similar stream function. 


W) 


1 

Re 2 


E 

y=i 


D i 


sin 6 Pj ( cos 0) 


(5.35) 


where Dj is a constant, and Pj ( cos $) is a first order associated Legendre 
polynomial. The far field dipole solution of (2.50) when combined with (2.66) 
gives 


G dipole (£ 0 ) 


1 sin 2 0 

4x f 


(5.36) 


With (5.36) we can calculate in (5.35) and find that it is 
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Re 2 

4ir 


(5.37) 


D i 


The D l coefficient we get for free, but the multipole coefficients (also known 
as tesserals) for D } - , where j > 2 require some extra effort in the form of 
Fourier integrals and the integral form of the Poisson equation. Because the first 
order associated Legendre polynomial is a complete orthogonal function, we may 
recast (5.35) by Fourier theory into the following integral 

7T 

A = Re 2 ? 1 ( ) Jf^j- f O ({,*)/>,'(«« f)df .(5.38) 

The self similar stream function G(£,0) seen in (5.38) is found by starting 
with an expression for vorticity by combining (2.18) with (2.7) and (2.4) giving 

V(v -5)-v 2 ^ = ^ • (5.39) 

To go beyond this point we follow Batchelor [1967] in recognizing that if the 
vorticity normal to the surface of the control volume containing our flow is equal 
to zero, then the following gauge (a Columb gauge) is appropriate 

V ’ S = 0 . (5.40) 

Equation (5.39) combined with (5.40) is now a Poisson equation which has the 
following integral form 

5 = — f 3 (*' > dV 

4tt J | 

control volume ‘ J 


(5.42) 
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where t is a position vector within the control volume. In the axisymmetric 
case, B has only one component B ^ which is related to ^(£,0) via (2.30) 
and (2.62). One then replaces the G(£,0) of (5.38) with (5.42) and proceeds to 
integrate. One can integrate the resultant equation by brute force replacing 
|7- af | with a “law of cosines” expression leading to relations involving ellip- 
tic integrals, or one can replace |if- af | with a Green’s function 


1 1 - * 


'* I 


00 l J l 

E E 

/=0 m=-l ^ 


ILz ml 

(/+m)« 


FT {cos P? (cos 0) e {m ^* ' ) 


(5.43) 


We must employ the axisymmetric form of vorticity and vector potential 


u5(af ) = - Uf sin <f> 1 e x + cos <j> e y (5.44) 

= — ~—r I - sin <(> e. + cos <f> e l (5.45) 

rsin 9 L "J 

where e x and e y are Cartesian unit vectors. Equations (5.43) - (5.45) are 
inserted into (5.42) and the real part is taken. The result is then inserted into 
(5.38) and the stream function and vorticity are converted to their self similar 
form via (2.66) - (2.67). The resultant equation is integrated employing the 
orthogonality relations for cos <f> and P}(cos 9) giving 
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The integration of (5.46) requires that either an analytic form of vorticity be 
known or a table of numerical values be provided for numerical integration. 
Equation (5.46) is the equation actually used in the software to find D t , 
where (5.35) could be used via (5.10) and (5.11) to define a Dirichlet boundary 
condition. The parameter Zoo is the £ value where the multipole boundary 
condition is imposed and presupposes that the vorticity for £ > £co * s 
sufficiently small as to be negligible. Only the first six terms of the series are 
used (1 < 6) with the belief that extra terms would be buried in the numerical 
error. 

It was later found that greater stability in the solution could be achieved if a 
Neumann rather than Dirichlet boundary condition was used for the stream func- 
tion at Zoo . This was brought about by restructuring the finite difference 
equation (5.26) for the n th node points in an n X n mesh as follows: 
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where 
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The partial derivatives in (5.47) are found analytically from (5.13) and (5.35). 
The vorticity boundary condition for was kept as a Dirichlet type with 

y(£ , 6 ) = 0 . Going to the Neumann boundary condition on stream 

function, Z enhanced the convergence rate and significantly reduced oscillation 
in the far field solution due to boundary condition mismatch. 

Equation (5.46) has a nice side benefit in that over all solution accuracy 
can be measured by comparing the Pj value calculated by (5.46), to the ana- 
lytic value given by (5.37). The D x value actually used in the boundary con- 
dition of (5.35) was the value given by (5.37). The numerical calculation of 
(5.46) involved using the latest iteration of IV(£,0) and integrating by 
Simpson’s rule in the f direction and in the 6 direction. The integrating 
formula used was the following, showing the f direction integral as an exam- 
ple: 

h 1 7 27 9 

/ f{y)dy = J /o + 4f/+ 2£+ ~ /„_2 + -g“ ( fn-l + /») + fn+\ 

(5.48) 

where 

U = h+h+h + • • • + /n-3 
E = / 2 + f\ + U + ' ' * + fn-4 

fn + 1 / (£oo) 
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/o = no) 

Coo = h{n+ 1) 

For (5.48) there are n node points between the two boundary node points at 0 
and n + 1, where h is the step size between each node point. The Q 
integration is the same as (5.48) except it is the upper limit instead of 
and a metric term is introduced to account for the theta stretching algorithm. 

4. The Angular Coordinate Mesh Stretcher 

In (5.14) we showed that the angular coordinate 6 used in the numerical 
computation is based on 6 and Reynolds number. Up until now the 8 coordi- 
nate has been left in general terms. The purpose of this coordinate is to act as a 
mesh stretcher so that computational points may be transferred from regions 
with low gradients and low truncation error to regions where the gradients and 
truncation error are high. A similar stretcher for f has not been developed 
because more points can be bunched near the origin by simply reducing the value 
of and bringing the multipole boundary condition in closer. The discovery 

of the most appropriate theta stretcher was achieved by a long process of trial 
and error. The ultimate version arose from the realization that the stretcher 
needed two control “knobs” to be effective. One knob controlled the angle at 

which no stretching occurred. This angle is referred to as 6 with greatest 
concentration of node points near the axis of the jet. The other knob was X, 
“rate-of-stretch” which could be set to zero where 5 = 


0 or to some 
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number causing a rearrangement of -points but not affecting the influence of 
0 . The resultant stretcher is 


p T 

f 1 

1 1 + X( «-*» - 1) | 0 - Xjt 

1 

1 


(5.49a) 


0 = -— In 
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(5.49b) 


where 

S = the stretched coordinate in radians 

X = rate-of-stretch knob 

0 = the physical space angle in radians 

i ] = the switch-over angle parameter 

0 = switch-over angle . 


To find 0 “switch-over angle” (rads) one uses (5.49b) with tj as the single 
parameter. 


One unfortunate aspect of (5.49a) is that while it is straightforward to go 
from 9 to S, it is difficult to go from 6 to 0 because (5.49a) is transcen- 
dental. To go from <5 to 9 one uses the following iterating formula 
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where 0 n is the present guess of theta and ^n+l the next guess. This for- 
mula is based on Newton-Raphson and is quite stable, converging in approxi- 
mately ten iterations. The first guess for starting 0 n is 6 . In the equations 
where 6 appears explicitly, (5.50) is used to calculate its value. 

Equation (5.49a) has several advantages as a stretcher: 5 = 0 for 

0 = 0 just as 6 = 7 T for 0 = n no matter what value of X or q is 

selected thereby ensuring that the boundary conditions of 0 = 0, ;r are not 
disturbed by the stretcher. While the function is valid only for the region of 
0<9<k, 0<6<ir it is singled valued and continuous in all derivatives. 

The 0 marks the angle where node points are being removed from the param- 
eter region of larger angle and added to the region of smaller angle. The 
appropriate value for t) is readily found by realizing that all the major action 

of the flow must occur in the region of 0 < 0 < 0 . The region defined by 

0 is a cone where the flow geometry (in State 3) is a mushroom shape that 
must fit snugly inside this cone (see Fig. 5-3A). This stretch algorithm has been 
referred to as a Japanese Fan coordinate stretcher. In Fig. 5-3B we can see an 
unstretched mesh (radial mesh not shown). As the Reynolds number is increased, 

coordinate mesh lines are transferred over from one side of 0 to another caus- 
ing greater density near the axis at 0 = 0. This is shown in Figs. 5-3C and 5- 
3D for increasing Re. This is analogous to a closing Japanese fan where the den- 
sity of fan pleats increases as the fan closes. One of the fortunate observations of 
this research was that the spreading angle of the State 3 round jet along with the 





Figure 5-3A-D. The Theta Stretcher as a Japanese Fan. 
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q parameter of (p, q) space is a linear function of Reynolds number (contrary 
to the Re 4 dependence flow, see Eqs. 4.12 and 4.14). This assumes that the 
numerical method has not become ill conditioned (stiff), at which point it 
becomes Reynolds number invariant. This makes deriving an empirical formula 
for t) = f(Re) quite simple. All that is required is two data points and know- 
ing at what Reynolds number a given mesh will go stiff necessitating a finer 
mesh. The X “rate-of-stretch” parameter was likewise assumed to be a linear 
function of Reynolds number and was described by a similar empirical formula 
based on two data points which was refined with higher Reynolds number. X 
was tuned by first establishing the correct i} by knowing the basic flow 
geometry and tuning X to higher values until the region outside the “switch- 
over cone” started to oscillate unnaturally, at which point X was backed down 
until the whole flow field was uniformly smooth. One aspect of the theta 
stretcher that should be understood by the reader is that it was not activated 
until the Reynolds number was greater than 10. For Re *— 10 the value of 
X was zero and the stretcher was inactive, resulting in a mesh of constant step 
size. 

5. The Matrix Solver 

We can now address the method of actually solving the finite difference 
equations of (5.25) and (5.26), which employs a special linear system solver. 
The linear system to be solved is the classic one of Q X = Y, where Q is a 
square sparse matrix, and X and Y are column vectors. The matrix Q con- 
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tains the finite difference equation, X is the solution vector, and Y is the 
right-hand side of the finite difference equation plus boundary condition terms. 
Since both .finite difference equations (5.25) and (5.26) are two-dimensional, 
second-order central difference equations, the molecule or stencil used to solve 
for one node point is made up of five points, see Fig. 5-3. The linear system that 
will solve a mesh based on this five-point molecule is the tri-block system which 
is analogous to the tri-diagonal system except that each of the elements along the 



Fig. 5-4. Control Difference Molecule 





diagonal are made up of matrices rather than scalar elements. The arrangement 
of this system is shown in (5.51) 



The subblocks of A,-, Bj, B? % Y, have the following form: 
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where k = 1 or 2, and a {{ , a ii±l , x,y, y,y, b- are scalar matrix elements. 
The matrix subblock Ay of (5.52) is a tri-diagonal matrix while Bj of (5.53) 
is a simple diagonal matrix. Using (5.26) as an example we can demonstrate how 
a finite difference equation is loaded in this linear system, which in this example, 
is centered on the node point at i ,j. For the j subblock 

= N - B 

a, ,• = C — 2(N + D ) 

= N +B 


bl = D +E 



For the j - 1 subblock, 


*;,• = d - e 

The computational domain is made up of n X n points (not including the boun- 
daries), see Fig. 5-1. One can see that for the case of i = 1 the Landau Squire 
boundary has to be taken into account due to i - 1, j component of the finite 
difference equation. This is done by setting the element for the Aj 

subblock equal to zero and taking the appropriate form of the Landau Squire 

solution multiplied by - (N - B ) and added to y,y . The multipole boun- 
dary is treated similarly while the 0 = 0,t boundaries result in the 

appropriate b £• elements being zeroed. Equation (5.51) is solved by getting it 
into an upper triangular form and solving by back substitution. The matrix in 
upper triangular form U X = R is the following: 



(5-54) 


where I is the n X n identity matrix. The equations relating (5.54) with (5.51) 
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are first the “boot strap equation” 


£/, = [(*)-' A,]"' (5.55) 

A, R, = Y, 

Then the “downward sweep equations”: 



(5.57) 

v, = 

(5.58) 

LjR, = Yj-B^Rf., 

(5.59) 

There is then the “backward sweep bootstrap” 

as 

II 

(5.60) 

with finally the “backward sweep equation” 
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(5.61) 


Equations (5.60) and (5.61) provide the solution. There is a problem with this 
method involving the storage of the Uj matrices. The storage requirements in 
the Q X = F system, not including zeros and the solution vector X are 
given by the equation 


# of elements = 6 n 2 - 4n 


(5.62) 


For the ILY = R system, not including zeros, l’s and the solution vector is 
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# of elements = n 3 (5.63) 

If one is using a mesh of 60 X 60 points where n = 60 then the storage 
required by Q X = Y is 21,360 while the storage required by U X=R 
is 216,000. This order of magnitude greater storage requirement is a consequence 
of the Uj submatrices being nonsparse because they are the result of matrix 
inversions. The order of magnitude greater storage requirement does not need to 
be accepted because the U X — R system is only an intermediate step and 
contains the same information as Q X= Y . The answer to this dilemma is 
to keep only the Uj_ x submatrix when calculating Lj in (5.57) and to throw 
Uj_ i away afterwards. When the downward sweep is complete the Uj 

matrices will have to be regenerated in order to provide the solution vector in 
(5.61). This requires the following additional equation 

Uj-l = ( B }., ) ( Aj - Uj 1 ) (5.64) 

This equation, called the “yo-yo sweep” equation is actually just a backward 
sweep equation that is a companion to (5.61). Using (5.64), Uj_ x is generated 
for (5.61) and the old Uj is thrown away. This technique of regenerating Uj 
by the yo-yo sweep is essentially trading off a factor of two increase in computer 
time for an order of magnitude increase in storage. There is a difficulty in (5.64) 
however, which is that this equation appears to amplify round-off error. 

The mechanism of this round-off error amplification is not fully understood 
and deserves further study. For the purposes of this problem a fix was found by 


c-0. 
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periodically “refreshing” the Uj matrix by using a stored version rather than 
regenerating it by (5.64). The frequency of refresh was found by taking a linear 
system with a known solution vector whose Q matrix was ill conditioned. The 
solution vector calculated was then closely observed and whenever the roundoff 
error became significant (i.e., observable) the U } - submatrix would be refreshed. 

The maximum matrix ever calculated by this technique was for an 80 X 80 
computation mesh which corresponds to a 6400 X 6400 element Q matrix 
(all elements including zeros). 

6. The Under Relaxation Method 

Next we examine the under relaxation method used to solve for vorticity and 
stream function. This was the price we had to pay for the simplicity of a two- 
step (vorticity solved then stream function) solution method. The under relaxa- 
tion used was 


= G. + T ( G* +1 - G„] (5.65) 


T = 1 

0 < Re < 5.0 

(5.66) 

T = 0.7 

5.0 < Re < 6.0 

(5.67) 

T = 0.1 

6.0 < Re 

(5.68) 


where 


G 


n 


the previous iteration stream function matrix; 
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£»+l — the latest calculated stream function; 

G n+1 — the under relaxed stream function; 

T = the under relaxation scalar whose values are shown in (5.66) to 
(5.68). 

The region of 5.0 < Re < 6.0 encloses the first transition where, due to the 
critical point bifurcation, it was necessary to under relax. For T = 1 (5.65) 

is made inactive. So for Re < 5.0 there is no under relaxation. Under relaxa- 
tion, when it was in force, would occur with every iteration of the two-step 
cycle. 

7. Solution Method Overview 

The software is cold started at a low Reynolds (typically Re = 0.1) and the 
solver is primed with the creeping solution for the nonlinear scale functions 


Y ( 6 ) = sin 0 

(5.69) 

Z(0) = sin 2 0 

(5.70) 


With these scale functions as a first iterate the next iteration is found for the vor- 
ticity scaling function Y i} - through (5.25) via the linear solver with Z,- • held 
fixed. With the updated vorticity the multipole constants are also updated via 
(5. -16). The vorticity solution is then used with (5.26) to solve for the stream 
function scaling function Z {] - by use of the linear solver. With the new stream 


/ 
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function Z, ; - the vorticity is recomputed as before through (5.25). The new 
vorticity is compared against the previous iterate. If the difference is sufficiently 
small between two successive iterations then the solution is considered con- 
verged. If not, then the cycle repeats until convergence or until the iteration 
limit is reached and the program aborted. In the case of convergence the values 

of vorticity, stream function, and multipole constants are written onto disk. The 

i 

Reynolds number is usually incremented by 1. The Landau-Squire constant is 
recomputed along with the mesh stretching parameters and under-relaxation con- 
stant. The previous solution is used as a priming first iterate and the cycle is 

i continued until convergence for the new Reynolds number. This process contin- 

ues from a creeping Reynolds number of 0.1, incrementing first by 0.9 and 
then by units of 1 until the Reynolds number is sufficiently high that acceptable 
accuracy with a given mesh is no longer attainable. Software validation for this 
system is relatively straightforward. The “nonlinear scale functions” (in this case 
a misnomer) for the linear solution are now in equations (5.69) and (5.70). If the 
convective terms in the software are disabled and a creeping form of the 
Landau-Squire solution is used as a boundary condition, then the resultant solu- 
tion will always be (5.69) and (5.70). Therefore it is a simple matter of com- 
parison to see if the software is error free. Likewise the far-field boundary condi- 
tion can be replaced by the Landau-Squire solution (which is also used in the 
near field) for a test of the validity of a solution which includes the convective 
terms. For this case the Landau-Squire solution would be the observed result for 
the entire computational domain. This makes error detection a matter of com- 
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parison with an analytic result. Because of these tests we have a high level of 
confidence in the validity of our method and software. 


8. Some Unsuccessful Approaches 

Finally, we would not want the reader to be left with the impression that 
the method described was found without going down some deadends, or without 
realizing that there are better approaches than the one used. A deadend 
approach that cost us a lot of time made use of a formulation in terms of the fol- 
lowing equations: 
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where rj = cos 0 . 


(5.72) 


Equations (5.71) and (5.72) represent a compact form of equations (5.1) and 
(5.2). Because of their symmetrical form it is tempting to think that an analytic 
solution can be found (though we have not found it). The deadend that we 
encountered involved trying to convert (5.71) and (5.72) into a single central 
difference equation while retaining the simple structure of the equations. We 
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found out that the multi-valued nature of t; = cos 6 near 0 — 0 and k plus 
the need for a 13 point finite difference molecule created unavoidable problems in 
the boundary conditions. There were two alternative approaches to the numeri- 
cal solution that were considered but not used. One of these two alternate 
schemes would have involved inserting (5.72) into (5.71) and expanding the 
result into a finite difference equation with stream function as the single depen- 
dent variable. This finite difference equation would then be solved by a sym- 
metric penta-diagonal block solver using Newtonian iteration. This approach has 
a number of advantages: one pass solver, no relaxation constants, single boun- 
dary conditions, etc. The main drawback is that the finite difference equation is 
so long that it is very difficult to avoid transcription errors without having access 
to an algebra program such as MACSYMA. The second scheme involved writing 
(5.71) and (5.72) together in a single linear system with the stream function and 
vorticity stacked in a single solution vector. This method had the advantages of 
having relatively straightforward difference equations, no relaxation constant and 
boundary conditions on a single variable (stream function). The disadvantage is 
that the linear algebra in the matrix solver is very messy requiring a customized 
algorithm for the unsymmetric diagonals. However, when we look back with the 
clear vision of hindsight, this second scheme might have been the best way to go. 
The customized solver would have been trouble. However, we too had to write a 
special solver for our tri-block system anyway, so the added complexity of the 
custom soiver would not have been that great. The customized single pass solver 
without relaxation would have certainly been faster. If the reader intends to 
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investigate a problem similar to ours it is recommended that the second scheme 
be investigated first before reproducing our actually used approach or resorting to 
the first scheme. 
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Chapter VI 

COMPUTATIONAL RESULTS 

1. Flows Computed By The Program MAVIN 

The numerical method described in Chapter Y was developed in a program 
named MAVIN. The source listings of all the software developed in this work, 
including MAVIN, are presented in App. C. MAVIN was run on the Control 
Data Corp. computer, CDC 7600 at NASA Ames Research center. MAVIN was 
used in two basic grid configurations. The first grid configuration was for a 30 
point by 30 point mesh (n = 30) with a far field boundary radius (£ c<) ) set at 
15. The results for the 30 X 30 cases displayed in App. B started with Re = 
0.1. Then Re = 1.0 was computed, thereafter the Reynolds number was 
increased by units of 1 until Re = 30 was reached. Plots for the 30 cases 
appear in App. B, Fig. B-l (Re = 0.1) up to Fig. B-30 (Re = 30). The second 
grid configuration used the finer mesh of 60 points by 60 points, with the same 
far field boundary radius at = 15. The 60 X 60 mesh led to more accu- 
rate results than the 30 X 30 mesh but suffered from the usual tradeoff of being 
much more expensive to run. As a consequence, only 10 Reynolds numbers 
were computed, specifically Re = 4, 6, 10, 15, 17, 20, 23, 25, 27 and 30. 
The results for these cases appear in App. B in Figs. B-31 to B-40. For each 
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Reynolds number MAV1N generated a matrix representing vorticity, and a 
matrix for stream function. MAVIN also generated six multipole constants used 
in the far field boundary condition, the far field boundary location the 

mesh size n, the Reynolds number, and the theta stretcher beta parameter. 

2. The Arch Spline Used in Data Reduction 

The data were outputted as ASCII format records, transferred to disk 
storage on a VAX computer and archived on magnetic tape. The discrete points 
represented in the n X n matrices of vorticity and stream function had to be 
transformed into a continuous polynomial format for subsequent processing. This 
was achieved by applying the following spline polynomial to the data 

/(£, 0) = £ 2 ( + CoO 2 + c 3 0 + c 4 ) + f ( c b (ft + c 6 0 2 + c 7 9 + c 8 ) 

+ Cq0^ + Cj q(P + Cj + Cj2 (6.1) 

where / (f,0) can represent either vorticity or stream function and the c k 
are constants found by fitting the spline technique polynomial to the MAVIN 
matrices. The spline used is represented in Figs. 6-1, 6-2 and 6-3, which show 
a hypothetical case of a 6 point by 6 point mesh (n = 6). The two dimensional 
spline function is a second-order polynomial in the f direction and a third- 
order polynomial in the angular (0 or 5) direction. Figure 6-1 shows how each 
spline passes through three points in the radial direction (for fixed 0). One can 
see how for Region 1 the spline is matched to the value of node 0, represented by 
D x (Dirichlet match). Nodes 0 and 7 are boundary condition node points for 
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this 6x6 mesh. The spline from node 0 is then matched at node 1 and node 
matched at node 1 and node 2. The region 2 spline has its values matched at 
nodes 1, 2, and 3 as represented by D 2 . 



Figure 6-1 Radial Spline. 
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Figure 6.3 Spline Domains. 


The radial splines are only matched in actual value, not in slope. One should 
note that the region of validity for splines 2-5 starts and ends halfway between 
node points, leaving much of the spline unused. 

Figure 6-2 represents how the splines are set up in the angular direction. It 
shows node points with constant angular stepsize which is consistent with the 
6 unstretched coordinate [9 is the physical angular coordinate). At node 0 
(boundary node) of Region 1 the spline is matched to the value of the node (Diri- 
chlet) and also to the value of its slope (Neumann). This is represented by 
( N/D)i . The spline of Region 1 is also Dirichlet matched through nodes 1 and 
2. The next spline for Region 2 is both Dirichlet and Neumann matched 
(N /D ) 2 to spline 1 at node 1. This process is continued along the circle of con- 
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stant f . This same process occurs at node 7 for spline 7 except in the opposite 
direction. The two systems of splines meet at region 4 where a keystone spline 
is matched to nodes 2 - 5. This is purely a Dirichlet match which is done only at 
the “keystone” for purposes of closing the system. We see that in all the splines 
only a part of the spline is actually used in computation. Figure 6-3 shows the 
two-dimensional spline regions with respect to the MAVIN node points. Some 
sample MAVIN node points are in parentheses. Those nodes which have values 
of 0 or 7 correspond to boundary nodes. Once a matrix of c k values for each 
spline region has been calculated, then manipulation of vorticity and stream 
function becomes quite straightforward. If one wishes to calculate a particular 
vorticity, W(f, 0), the f coordinate is divided by h stepsize to find the i 
index and 6 is converted to 8 . The value of 8 is then divided by the k 

stepsize to yield the j index. With the i,j indices known, the appropriate 
spline region is selected and the c k values are pulled from its storage matrix 
and inserted into (6.1). Coordinates f and 0 are then inserted into (6.1) to 
yield the vorticity. This procedure can yield a higher resolution of the computa- 
tional domain. It can also provide a simple analytic equation in (6.1) for later 
manipulation in solving the particle path equations and determining (p , q ) 
values and critical point locations. 

3. Plots Of Computational Results 

The software that performs this spline operation is called INVORT and is 
shown in App. C. A higher resolution matrix was generated from these splines 
and digested by the contour plot programs CONTOUR and CONVORT and the 
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particle path plot program PLOT which produced the entrainment diagrams. 
These plots are displayed in App. B in self similar coordinates. Using the results 
shown in Figs. B-31 and B-32 we have produced Fig. 6-4 in which the particle 
path plot has been superimposed on the entrainment diagram. The momentum 
source is always at location (0.0, 0.0) on the plots shown in this work. The 
flow is always going from left to right along the x axis. The x axis is the axis 
of symmetry with the plots shown corresponding to the upper half of a plane 
passing through the axis of symmetry of the flow. The unshown lower half would 
simply be a mirror image. In Fig. 6-4 we see that the two stream function con- 
tour plots are basically dipoles with a bias slightly to the right of the origin. The 
following coordinates are used: 


x = 



( 6 . 2 ) 



(6.3) 


e = 


n/z 2 + j / 2 


(6.4) 


where x and y are physical Cartesian coordinates. Coordinates x and y are 
self similar coordinates and are the type used in the axis labels of App. B. We 
see that the center of the stream function contours for Re = 4 is at about (x, y) 
= (0.2, 2.0) and for Re = 6 at about (0.5, 1.8). This center corresponds to the 
peak value of stream function (note that the stream function is zero at the origin 
and at infinity). The vorticity contour plots of Fig. 6-4 are very similar to each 
other except for a slight bias to the right at Re = 6 as compared to Re = 4. 
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Fig. 6-4: Computed solutions for the round jet at a) Re = 4.0, and 

b) Re = 6.0. Quantities displayed are self-similar stream 
function, vorticity, and particle paths. 
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The peak value of vorticity is always at the origin since vorticity has a singular- 
ity proportional to l/£ 2 . The entrainment diagrams shown in this work are all 
based on the particle path equations (4.1) and (4.2). The particle path equations 
are used to produce the array of isoclines (all normalized to the same length) 
shown in the entrainment diagrams of Fig. 6-4 and App. B. In contrast to the 
vorticity and stream function just discussed, we see in Fig. 6-4 that the two 
entrainment diagrams at Re = 4 and 6 are significantly different. In the Re = 
4 case, one critical point is on the x axis at about ( x e , y e ) = (1.3, 0). In the 
Re = 6 case, there are two critical points. One is on the x axis at 
(x e , y e ) = (2.0, 0) and the other is off axis at (x c , y e ) = (1.6, 0.8). At this 

stage a key aspect of the work becomes apparent. For different Reynolds number 
the entrainment diagram topology can be completely different while the vorticity 
and stream function plots can be almost identical. In this respect, entrainment 
diagrams are a much more effective means of displaying the structure of the flow 
than are plots of stream function or vorticity. 

4. Transition In The Numerical Solution Of The Navier Stokes Equa- 
tions 

.As previously discussed in Chapter 3 the bifurcation from one topology to 
another represents a form of transition. Based on the results for the p , q 
plot discussed (in Ch. 3) we know that Re = 4 must be in State 1 since it has 
one on-axis critical point (which has to be a stable node). A close examination of 
Fig. B-41 shows that the entrainment diagram vectors are consistent with a 
stable node for Re = 4. We can bracket the first transition Reynolds number 



- 94- 


Re! by studying Fig. B-5 for Re = 5 where again we see only one critical 
point thereby concluding that Re = 5 is also in State 1. We now know that 5 
< Re j < 6. A naive way of refining Re x would be to run a series of compu- 
tations for Re = 5 to 6 looking for critical point bifurcation. This method is 
ineffective because the region near the axis is one of very high gradients and con- 
sequently high truncation error. A better way is to take values of the off-axis 
critical point (x e , y e ) for different Re and extrapolate back to Re v The 
parameters of the off-axis critical point are acquired through the program named 
RAMMER. HAMMER evaluates the particle path equations (4.1) and (4.2) at 
the off-axis critical point where the equations are equal to zero. The spline poly- 
nomial of (6.1) is inserted into these particle path equations yielding the following 
critical point location formulas: 


i_ r c _M + + c 8 1 

2 c l(^) + c 2(0f) + c 3 (0«) + c 4 


-f sin 0 c ~ 
2 c 

+ <7^ + C n 



c lfe + C b€c + c fl j 

+ 2 

c 2 fc + <*6^ + <70 


These transcendental formulas are iteratively convergent provided that the 
critical point location (f e , $ e ) is within the spline domain of the c k con- 
stants used in (6.5) and (6.6). HAMMER works interactively with its user by 
first receiving a guess of the critical point location. The first guess is used to 
select a spline domain and a second guess is found by (6.5) and (6.6). If the 
second guess is within the original domain then the critical point is captured. If 
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not, a new domain is selected based on the second guess and the process is 
repeated. If the critical point is on the corner between different spline domains, 
the HAMMER software can go into an endless hysteresis loop so the user is con- 
tinually given updated information on domain shifts. Breaking this hysteresis 
loop requires only a simple software adjustment in the spline domain selector to 
bias HAMMER toward a particular choice. Once the critical pont is captured, its 
coordinates and q value are calculated based on the current spline. In Table 6-1 
the critical point parameters based on the 30 X 30 mesh (n = 30) results of 
App. B are shown. In Table 6-2 a similar table is shown for the 60 X 60 mesh 
(n = 60) results for comparison of numerical effects. From Table 6-1 we can 
extrapolate the first transition Reynolds number Re t by fitting a parabola 

through Reynolds numbers 6, 7, and 8 with q as the independent variable. 
This gives 


Re = -2.22q 2 + 5.28q + 5.35 (6.7) 

We know from p , q theory that the first transition occurs at q = 0. We can 
extract Re t from (6.7) by letting q = 0 giving Re! = 5.35. This result 
reflects a particular set of data and numerical method. The naive approach men- 
tioned earlier of fine scanning between Re = 5 through 6 yielded Re!=5.1. 
This value can be regarded as a lower limit. Of all the different data bases exam- 
ined in the development of MAVIN, the highest value ever observed was 
Re ! = 5.9. Thus a conservative evaluation of the first transition Reynolds 
number is Re t = 5.5 ± 0.4 with 5.4 being the value most consistent with 
the App. B data. Before leaving Fig. 6-4 we should note that the off-axis critical 
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Table 6-1 

CRITICAL POINT PARAMETERS FROM APPENDIX B, N = 30 
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Table 6-2 

CRITICAL POINT PARAMETERS FROM APPENDIX B, N — 60 



Off-Axis Critical Points 

On- Axis Critical Points 

Re 

*c 

Vc 

o c 

q 



N/A 




1.35 

6 

1.61 

0.80 

26.5 

0.17 

2.00 

10 

1.98 

1.29 

33.2 

1.27 

3.25 

15 

3.06 

1.48 

25.8 

2.80 

4.75 

17 

3.48 

1.52 

23.7 

3.13 

5.30 

20 

4.03 

1.58 

21.3 

4.02 

6.05 

23 

4.49 

1.62 

19.8 

4.25 

6.65 

25 

4.77 

1.64 

19.0 

5.14 

7.05 

27 

5.00 

1.65 

18.3 

4.94 

7.40 

30 

5.30 

1.67 

17.5 

6.01 

7.80 
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point location for Re = 6,0 was at (1.61, 0.80), while the stream function 
contour center was at (0.5, 1.8). This emphasizes that there really is no simple 
correlation between the off-axis critical point and the stream function maximum 
point. The off-axis critical point for Re = 6 has a q value of 0.17. This 
value is consistent with a stable node (0 < q < 9/16) thereby showing that 
Re = 6 is in State 2. The flow will undergo the second transition into State 3 
when q attains and surpasses a value of q = 9/16 = 0.56. We see from 
Table 6-1 that this second transition Reynolds number has the following limits: 
7 < Re 2 < 8 . A high resolution scan was performed through this region yield- 
ing the data shown in Table 6-3. A parabola was fitted through Re = 7.4, 7.5 
and 7.6 with q as the independent variable giving 

Re = 11.90q 2 - 10.12q +9.46 . (6.8) 

With q replaced by 9/16 we find Re-j = 7.54. In comparing Table 6-1 

Table 6-3 

HIGH RESOLUTION SCAN THROUGH Re 2 


Re 

*c 

Vc 

8c 

q 

7.4 

1.62 

1.09 

33.9 

0.51 

7.5 

1.63 

1.10 

34.0 

0.55 

7.6 

1.64 

1.11 

34.2 

0.58 

7.7 

1.84 

1.12 

34.4 

0.62 
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with Table 6-2 we find a q variation of 0.04. Assuming a q deviation of 
0.06 for Re 2 we can determine a standard deviation from Table 6-3 for 
Re 2 of 0.16. This gives our final result as Re 2 = 7.54 ± 0.16. 

5. The State Three Solution 

For Re > Re 2 the flow goes into State 3 where the off-axis critical point 
becomes a stable focus. In Fig. B-8 we observe this stable focus just beginning to 
appear in the entrainment diagram. In Fig. 6-5 (Re = 15) a particle path tra- 
jectory has been drawn over the entrainment diagram. For this case, the stable 
focus is well established and we see the characteristic mushroom shape of the 
round jet. The critical point on the axis is a saddle point as seen in Fig. B-41 
for the case of Re = 10. In Fig. 6-5 the center of the stream function contour is 
at about (x, y) = (2.3, 2.25) while the off-axis critical point is at (3.06, 1.48). 
This further emphasizes the dissimilarity between the entrainment diagram and 
stream function contour plot. It is important to note that on the vorticity con- 
tour plot there is no local concentration of vorticity at the critical point location 
(3.06, 1.48). In fact, this region in the vorticity field is smooth and decreasing 
monotonically with increasing radius with no indication that a stable focus has 
formed. We can observe in both Figs. 6-5 and 6-6 that the vorticity and stream 
function contours are shifting more to the right as Reynolds number increases. 
The effect of increased Reynolds number on the particle paths as seen in Fig. 6-6 
is that the “stem” of the mushroom shape is lengthening and the intensity of the 
stable focus is increasing. In Fig. 6-7 the off-axis critical point angle 0 C vs 
Reynolds number can be observed. In the State 2 region of Re j < Re < Re 2 
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the flow is very dynamic with a rapidly increasing critical point angle with slowly 
increasing Reynolds number. The maximum 6 e of 34.50 for Re = 8.5 is 
reached shortly after the second transition into State 3. After this peak 6 e the 
critical point angle decreases slowly. Intuitively one would anticipate this angle 
getting smaller for larger Reynolds number since one might expect the mushroom 
stem to grow faster than the diameter of the mushroom. 

6. Numerical Instability 

In Fig. 6-7 we run into numerical difficulties at about Re = 15 where the 
30 X 30 mesh and 60 X 60 mesh start to yield different results. This is further 
shown in Fig. 6-8 where q is plotted vs Reynolds number. Figure 6-8 
represents one of the biggest surprises in the study. For q to have this linear 
behavior is totally unexpected. Low Reynolds number theory predicts an Re 4 
dependence for q. We see in Fig. 6-8 that Re = 15 is a point of divergence 
between the two meshes. 

If App. B is examined for the 30 X 30 mesh cases, one can see the manifes- 
tation of numerical error in the formation of vorticity clumping, oscillations in 
the stream function and dimpling at the saddle point in the entrainment 
diagram. These phenomena are completely numerical in origin. This can be seen 
in Fig. 6-9 where two vorticity and stream function contours are compared with 
the only difference being the mesh resolution. The vorticity clumps for the n = 
30 case are centered around every other node point in the radial direction. 




OFF-AXIS CRITICAL POINT ANGLE 
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Fig. 8-7. Off Axis Critical Point Angle vs Reynolds Number. 
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Flg. 6-8 


q vs Reynolds Number 
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This clumping completely disappears in the n = 60 case. For the stream func- 
tion we can see how the contours shift to the right and become free of oscillation 
for the finer mesh. The effect of numerics is even more pronounced in Fig. 6-10. 
Here we compare the dipole coefficient D x (which is the first coefficient in the 
multipole far field boundary) for n = 30, 60 and the exact result, versus Rey- 
nolds number. Since we have an exact result in (5.37), we know at what Rey- 
nolds number the solution is starting to degrade. It can be seen that the curves 
are starting to separate at Re = 12. The n = 30 curve changes shape alto- 
gether between Re = 21 and 22. Examination of Fig. B-21 which has the con- 
tour plots of vorticity and vorticity times I ; 3 (which is the kernel function used 
in solving D x ) shows that clumping is just beginning at Re = 21. The “dog 
leg” of n = 30 in Fig. 6-10 is never observed for n = 60, nor has vorticity 
clumping or the other overt symptoms been observed. However it is probable 
that if we had carried the n = 60 calculations much beyond Re = 30 these 
features w'ould have been observed. A question that immediately arises is how 
accurately does an error between the computational and exact value of D x 
compare with error in the vorticity distribution. This point is emphasized by the 
plots of vorticity (W) and vorticity times f 3 , (VF*f 3 ) shown in App. B. 

From the App. B plots we see that W>f 3 starts to clump before W, which 
would lead one to believe that D x may be an overly sensitive measure of accu- 
racy since Wf 3 is the 



ORIGINAL PAiiL * * 

OF POOR QUALFTY 


STREAM FUNCTION CONTOUR PLOT FOR Re-30.0 

30 X 30 MESH FOR (.=15,0 


VORTICITY CONTOUR PLOT FOR Re-30.0 

30 X 30 MESH FOR t =15.0 


STREAM FUNCTION CONTOUR PLOT FOR Re-30.0 

60 X 60 MESH FOR £..=150 


VORTICITY CONTOUR PLOT FOR Re=30.0 

60 X 60 MESH FOR £ =15.0 







Computed Solutions for the round jet at Re = 30 showing: 
a) vorticity clumping due to numerical instability on a 30 
X 30 mesh; b) smooth solution on a 60 X 60 mesh. 
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Fig. 6-10. Dj Tesseral Compared vs Reynolds Number. 
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kernel in the D l calculation. We can determine this degree of sensitivity from 
(5.37) for the dipole term which is 

Coo jt 

D x = f f W(£ ' , 9 ' )t' 3 sin 2 0 ' do' d£' (6.9) 

o o 

We will assume that a vorticity error AW will produce a dipole error of 
AZ)j. We can write down vorticity as 

W[Z, 0) = W(£, 0) + AW (6.10) 

where W(£, 9) is the computed vorticity and W(f, 9) is the actual. We may 
expand both W(f, 9) and AW into their multipole constituents. However, 
by orthogonality we know that only the dipole components will contribute to 
(6.9). We can therefore write down vorticity as 

W(f, 9) = sin 9 [ Ri{Z) + Ai?j j + ••• (6.11) 

where R x is the true vorticity radial dipole component and A R x the radial 

dipole error component with higher-order terms dropped. Equation (6.11) is 
inserted into (6.9) and the angular integration is performed giving 

*oo 

A + AD, = ££ f [«,(«') + AT?,] e' 3 #' (0.1-21 

0 

If we assume that A/?j is a relative error such that £R l /R l is a constant 
and not a function of 
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D. + AA = ££[l + ^.] / *,(«'>«' 3 


(6.13) 


Since 


we conclude that 


D, = 


Re 2 


/ *,«')«' 3 «' 


(6.14) 


AZ), 

"a" 


A/?J 

"rT 


(6.15) 


However, if we assume that A/?j is an absolute error which is a a constant 
and not a function of £ then by combining (6.12) with (6.14) and integrating 
with A/?j outside the integral and dividing by (5.37) we find that 


A A AR£ to 4 tt 

D x 3 


(6.16) 


Note that the observed numeric error in vorticity appears to be somewhere 
between the two cases of absolute and relative error. The method used in 
MAVIN calculates a nonlinear scale function and then multiplies that result by 
the analytic representation of the radial component of the creeping flow. Since 
this analytic representation is dominated by a Gaussian function, we know that 
the ultimate vorticity function’s error diminishes with distance and is ultimately 
very small. 

There is round-off error and truncation error in MAVIN and in the quadra- 
ture used to calculate D x . However since the CDC 7600 uses a 60 bit word and 
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double precision was used in many of MAVIN’s computations we know that 
round-off error is insignificant. The chief culprit for error in this computation is 
truncation error. What we can conclude is that (6.15) represents a best quality 
factor expected in the vorticity field. Specifically if 

AD, 

1-— = 0.80 

then we must assume that the overall quality of the vorticity is 0.80 or worse. 
The main advantage of this quality parameter is that it tells us when to no 
longer trust a calculated result from MAVIN. We can employ this in examining 
the parameter 


Numeric D l 
Analytic 

in Tables 6.4 and 6.5. Our own rule of thumb is that quality of less than 0.80 
is not acceptable. Based on this the 30 X 30 mesh results are acceptable up to 
R e = 18 but not beyond. The 60 X 60 mesh results are acceptable up to 
around Re = 23. There is, however, a paradox arising from this analysis. For 
the mesh of 30 X 30 in the Re = 21 case we observe some vorticity clumping 
occurring in the \V*£ plot of Fig. B-21. The quality factor for B-21 is 0.73. 
However, in Fig. B-40 for the 60 X 60 case of Re = 30 we observe no clump- 
ing or any undesirable numeric effects, even though the quality is an abysmal 
0.67. A possible solution to this paradox can be offered through a line of specula- 
tion. The multipole coefficients or tesserals of Table 6-4 are plotted in Fig. 6-11. 
We observe with some surprise that all the coefficients come to almost a single 
point near Re = 8. 


Ill - 


Table 6*4 

MULTIPOLE COEFFICIENTS FOR A 30 X 30 MESH AT = 15 


Numeric Results 

_ Numeric D, 

Re — — - — ; — — 

Analytic 

D\ D% D+ Cj D$ 

0.1 1.00 0.796E-3 0.135E-6 -0.208E-6 0.797E-7 -0.788E-6 0.152E-5 

1 1.00 0.796E-1 0.130E-2 0.244E-5 0.543E-5 -0.776E-4 0.137E-3 

2 1.00 0.318E0 0.208E-1 0.140E-2 0.948E-4 -0.297E-3 0.379E-3 

3 1.00 0.716E0 0.105E0 0.165E-1 0.262E-2 -0.255E-3 0.295E-3 

4 1.00 0.127E1 0.327E0 0.912E-1 0.283E-1 0.651E-2 0.119E-2 

5 1.00 0.199E1 0.785E0 0.336E0 0.150E0 0.669E-1 0.274E-1 

6 1.00 0.286E1 0.158E1 0.949EO 0.598E0 0.385E0 0.245E0 

7 1.00 0.388E1 0.280E1 0.220E1 0.182E1 0.154E1 0.132E1 

8 0.99 0.502E1 0.451E1 0.442E1 0.456E1 0.484E1 0.524E1 

9 0.97 0.827E1 0.676E1 0.798E1 0.990E1 0.127E2 0.166E2 

10 0.96 0.762E1 0.958E1 0.132E2 0.192E2 0.288E2 0.443E2 

11 0.97 0.939E1 0.138E2 0.225E2 0.385E2 0.682E2 0.123E3 

12 0.96 0.110E2 0.182E2 0.334E2 0.644E2 0.129E3 0.263E3 

13 0.94 0.127E2 0.232E2 0.470E2 0.101E3 0.223E3 0.505E3 

14 0.91 0.143E2 0.286E2 0.633E2' 0.148E3 0.358E3 0.888E3 

15 0.89 0.160E2 0.343E2 0.818E2 0.206E3 0.539E3 0.144E4 
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Table 6-4 (Cont.) 


Numeric Results 

_ Numeric D 1 

Re — — ■ — . - 

Analytic D\ 

D l D% D, Dj Di Df 

16 0.86 0.176E2 0.404E2 0.103E2 0.277E3 0.773E3 0.221E4 

17 0.83 0.192E2 0.467E2 0.126E3 0.359E3 0.106E4 0.322E4 

18 0.80 0.207E2 0.529E2 0.150E3 0.450E3 0.140E4 0.446E4 

19 0.77 0.221E2 0.591E2 0.175E3 0.546E3 0.177E4 0.591E4 

20 0.74 0.236E2 0.656E2 0.202E3 0.65763 0.222E4 0.77 1E4- 

21 0.73 0.256E2 0.747E2 0.241E3 0.827E3 0.294E4 0.107E5 

22 0.70 0.269E2 0.808E2 0.270E3 0.949E3 0.347E4 0.131E5 

23 0.65 0.274E2 0.835E2 0.282E3 0.101E4 0.373E4 0.142E5 

24 0.61 0.278E2 0.853E2 0.291E3 0.105E4 0.392E4 0.151E5 

25 0.56 0.280E2 0.867E2 0.298E3 0.108E4 0.408E4 0.158E5 

26 0.52 0.282E2 0.879E2 0.304E3 0.11 1E4 0.422E4 0.165E5 

27 0.49 0.283E2 0.890E2 0.310E3 0.114E4 0.436E4 0.171E5 

28 0.46 0.285E2 0.900E2 0.31 5E3 0.117E4 0.450E4 0.178E5 

29 0.43 0.286E2 0.911E2 0.321E3 0.120E4 0.465E4 0.185E5 

30 0.40 0.288E2 0.926E2 0.329E3 0.124E4 0.485E4 0.195E5 
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Tftble fl-5 

MULTIPOLE COEFFICIENTS FOR A 60 X 60 MESH AT — 15 

Numeric Results 

Numeric D x 
Re 

Analytic D 1 

Dx D 2 D s D 4 Di D 6 

4 1.0 0.127E1 0.328E0 0.914E-1 0.265E-1 0.776E-2 0.196E-2 

6 1.0 0.286E1 0.159E1 0.956E0 0.604E0 0.392E0 0.258E0 

10 0.99 0.785E1 0.101E2 0.142E2 0.210E2 0.321E2 0.502E2 

15 0.96 0.172E2 0.390E2 0.977E2 0.259E3 0.71 1E3 0.200E4 

17 0.93 0.215E2 0.563E2 0.164E3 0.504E3 0.161E4 0.526E4 

20 0.88 0.280E2 0.867E2 0.298E3 0.109E4 0.412E4 0.160E5 

23 0.82 0.343E2 0.120E3 0.468E3 0.193E4 0.827E4 0.363E5 

25 0.77 0.384E2 0.144E3 0.598E3 0.264E4 0.121E5 0.565E5 

27 0.73 0.422E2 0.167E3 0.735E3 0.342E4 0.165E5 0.820E5 

30 0.67 0.478E2 0.203E3 0.955E3 0.476E4 0.246E5 0.131E6 














TESSERAL STRENGTH 



Fig. 6-11. Tesseral Strength vs Reynolds Number 
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Since all of these coefficients are sensitive to numerical error (compare D z in 
Tables 6-4 and 6-5), Fig. 6-11 is useful only in a qualitative sense. Our first 
point of speculation is that the Reynolds number for which the Dj coefficients 
become equal is Re 2 , the second transition Reynolds number. This seems 
intuitively reasonable since the flow at Re 2 consists of an off axis star point 
which suggests a relatively simple flow topology. Since all Dj coefficients 
must be equal to zero at Re = 0 we bring about our second line of speculation 
based on Fig. 6-11 that the Dj coefficients have the following analytic form 


D: 


1 RC 2) 

Re 

1 4ir ) 

Re 2 


(6.17) 


If (6.17) is inserted into the multipole stream function of (5.35) and the Cauchy 
ratio test is performed for 0 = jt/2 one finds that convergence occurs only if 


Ifl > \0 ■ (6.18| 

With tog defined as the lower limit of the multipole expansion we find 
that for = 15 the Reynolds number must be less than or equal to 29.22 
for convergence of the infinite multipole series. Clearly for a six term finite series 
approximation of the multipole infinite series, the Reynolds number would have 
to be much less than 29 for there to be good accuracy. What is suggested by 
this line of speculation is that the low quality of Re = 30 for 60 X 60 is not 
so much a reflection on the numerics but rather on the accuracy of the truncated 
multipole expansion when applied at = 15 for this Reynolds number. 
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7. The Computer Animation Of Axisymmetric Jets 

Perhaps the most interesting aspect in the data analysis of this study was 
the creation of an eleven minute computer animation. This movie was comprised 
of six parts: the first four parts are Navier Stokes solutions based on a 60 X 60 
mesh of the round jet at Re = 4, 6, 20 and 30 (States 1, 2, and 3). The last 
two parts are Stokes solutions of the vortex ring and ramp jet. Animation of the 
Stokes solutions was straightforward since their particle path equations are ana- 
lytic. Equations (4.3) and (4.4) were marched in time using a fourth order 
Runge-Kutta integrator to produce a time line for the vortex ring. The software, 
which is shown in App. C started the time lines as straight lines in physical coor- 
dinates. In Fig. 6-12 the momentum source is shown as a plus and the critical 
points are asterisks. The p , q plot is in the upper right-hand corner with the 
flow state and Reynolds number shown below. The initially straight time lines 
are drawn instantaneously at Re = 100 since the equations are singular at t = 
0 where Re is infinite. As seen in Fig. 6-12 the time line immediately starts rol- 
ling up around the stable focus forming the expected mushroom shape. In Fig. 
6-13a, the Reynolds number decreases to below Re 2 and the flow goes into 
State 2 where the rolling stops. The off-axis critical point then merges into the 
on-axis saddle and the flow topology changes to State 1. Figure 6-13b shows the 
entire past flow history from State 3 to the State 1 flow depicted which is 
effectively a dead flow. The p , q plot evolves as the time lines and critical 
points go through their motions, all of which are based on the equations of 
Chapter 4. The ramp jet involves the integration of (4.7) and (4.8). Unlike the 
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vortex ring, States 1 and 2 are not masked by the flow history as seen in Figs. 
6-14 through 6-16. Like the vortex ring the flow starts as two straight time lines 
in physical space. However, for the ramp jet the initial Reynolds number is 
small. As the Reynolds number increases the flow bifurcates from State 1 to 2 
and then to 3 with the corresponding transformation of the off-axis node to a 
stable focus. Figures B-42 through B-47 show the time line evolution for the 
three Stokes solutions of the vortex ring, round jet, and ramp jet. These plots 
(in physical space) were the prototypes for the movie software and show some 
details in the fluid strain that are less apparent in the movie. The source is at 
coordinate (0, 0) and is blowing to the right. The critical points are not shown 
in Figs. B-42 through B-47. The crosses shown in Fig. B-42 are individual points 
being marched by the Runge-Kutta routine. At Re = 69.9 the crosses are so 
closely packed that they appear as a single line. At Re = 13.19 the individual 
crosses can be seen in the front of the mushroom (the front being the right-most 
part near the axis of symmetry), which is an area of very high strain. The coat- 
tails of the flow are still very dense but are (or were) being sucked into the volute 
revolving around the stable focus. It should be emphasized that a vortex ring 
stops rolling up after transition from State 3 to 2 and becomes a dead flow as 
Re vanishes. However, even though the flow is dead the mushroom remains as a 
consequence of the flow’s history. Figure B-43 is the same flow as B-42 but the 
time line is further away from the source. The consequence of this is that the 
time line never forms a mushroom because transition occurred before this line 
could form a volute. Figure B-44 shows a Stokes round jet for Re = 2.0 which 






Fig. 6-16 Stokes Ramp Jet. State 3 
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is in Stills 1* W s sss the characteristic sniooth parabolic shape that the State 1 
flow has for all time. Figure B-45 is a Stokes round jet for Re = 8 which is in 
State 2. A characteristic corner is forming in the time line. State 2 flows never 
roll up to form a mushroom shape. Notice that the front has a high degree of 
strain unlike the coattails. Figure B-46 shows a Stokes round jet for Re = 20 
which corresponds to State 3. This flow is rolling up and is forming the charac- 
teristic mushroom shape. The three round jet solutions are displayed in dimen- 
sional units because the round jet has no length scale. The fluid simulated in 
the round jet is olive oil. Figure B-47 displays the Stokes ramp jet. The ramp 
jet, like the vortex ring, goes through all three states but its flow history does not 
mask its present state. Unfortunately (as shown in the movie), the ramp jet goes 
through States 1 and 2 so quickly that their features are almost unobservable 
(Re ! = 3.7, Re 2 = 5.8 for the ramp jet). The features of State 2 are distinc- 
tive only in the round jet where the Reynolds number is not a function of time. 

All of these Stokes flows have an unrealistic aspect in that their volutes do 
not translate downstream but loiter around the momentum source. This is a 
consequence of the elimination of the convective term in the creeping approxima- 
tion. When one sees the Stokes flow movie this aspect of no translation gives the 
flow an unnatural appearance. However, as mentioned before, four Navier Stokes 
flows for the round jet based on a 60 X 60 mesh were also animated. Animating 
these flows was much more difficult because discrete numeric data had to be 
integrated rather than analytic functions and the flows were convecting making a 
correct viewing frame more challenging. The particle path equations (4.1) and 
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(4.2) were used with a stream function spline based on (6.1). This particle path 
representation was valid only out to ^ <X) . Beyond the multipole solution 
was used via (4.1) and (4.2). The equations were integrated using the same 
Runge-Kutta method as used in the Stokes solutions. The Navier Stokes simula- 
tion uses fluid parameters consistent with olive oil (glycerine was represented in 
the Stokes flow). Only the upper half of the plane is calculated and then 
reflected over the axis of symmetry for the lower half in the graphics. The inner 
time line (closest to the momentum source) is made up of 700 points with 600 
points bunched in the 1/20 part of the time line closest to the axis of symmetry. 
The outer time line is made up of 300 points with 200 points bunched in the 1/8 
part of the time line closest to the axis of symmetry. This large amount of point 
concentration near the axis of symmetry is required by the high rate of strain 
that occurs in the front of the round jet. The movie framing rate is at 24 
frames per second. One can see that the amount of data generated is enormous. 
An 11 minute movie has 15,840 frames with each frame having 1,000 points 
and each point described by two single precision words. This equates to 31.7 
megabytes just for raw data storage. This data in turn, has to be converted to a 
graphics format that increases storage requirements by an even larger amount 
(over an order of magnitude). The storage requirement exceeded available com- 
puter storage necessitating that the movie be made in roughly two minute seg- 
ments and spliced together. The actual filming was done by first generating the 
graphics on disk and then transferring to magnetic tape. The magnetic tape was 
read by a Dicomed machine which built up each image on a high resolution black 
and white CRT. The visual image on the CRT was passed through a color wheel 
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into an animator’s movie camera. Each frame of film was triple exposed with 
different images passing through different filters of the color wheel. The colors 
used were red for the critical points, white for the source, and yellow and green 
for the time lines. Figures 6-17 through 6-20 show frames from the Navier Stokes 
movies. 

States 1 and 2 are similar in appearance to their Stokes counterparts. Fig- 
ure 6-17 shows the smooth parabolic shape of State 1 for Re = 4 with the time 
lines not touching each other. Figure 6-18 represents a State 2 flow for Re = 6. 
Near the axis the time lines have merged with each other and changed color 
where they overlay. The saddle and off-axis nodes are on the time lines. Notice 
that the State 2 kink is not at the off-axis node. As time progresses the kink 
gets closer to the off-axis node and presumably with infinite time they are super- 
imposed. Being in State 2 (like State 1) this flow will never roll up and form a 
mushroom shape. Because the Reynolds number regime in which State 2 occurs 
is so narrow (between Re = 5.5 and 7.5) an experimentalist would have to be 
specifically looking for State 2 in a very carefully controlled experiment. The 
discovery of State 2 and the creation of the State 2 animation was one of the 
major dividends of this study. Two State 3 flows are shown in Figs. 6-19 and 
6-20. Figure 6-20 is by far the most dramatic flow simulated. This Re = 30 
flow looks incredibly natural and has been mistaken by some of the movie’s audi- 
ence as experimental data. In fact, the movie looks better than it should judging 
from its dipole coefficient quality factor which is a paltry 0.67. Numerically 
speaking the Fig. 6-19 could serve as a basis for experimental verification unlike 
the Re = 30 which is less accurate. In both flows we notice the long natural 
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dtems of the mushroom shape. As time progresses these stems are sucked up into 
the volutes which are translating away from the source. The saddle point always 
has the time lines bunching about the round jet. In round jets the “bunch point” 
where time lines accumulate the saddle point is one and the same. This is not 
the caste in the nonautonomous flows like the vortex ring and ramp jet where the 
bunch point is separate though obviously associated with the saddle point on the 
flow axis. In the Re = 30 flow if one looks carefully near the saddle point just 
before the movie ends, one can see a slight dimple and the early effects of spline 
kinks due to the high rate of strain. All of the Navier Stokes movies exhibited a 
remarkable naturalness. 
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Fig. 6.17. Navier Stokes, Re = 4, State 1. 
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Fig. 6.18. Navier Stokes, Re = 6, State 2. 
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Fig. 6.20. Navier Stokes, Re = 30, State 3. 
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Chapter VII 
CONCLUSIONS 

Two basic sets of conclusions can be drawn from this work. One set is 
drawn from the numerical methods employed. The other set is drawn from the 
physics of the problem under study. 

1. On The Numerical Method 

The two-step numerical method involving the solution first of the stream 
function with vorticity fixed and then of the Poisson equation for the vorticity is 
widely used and works quite well. However, it is the author’s opinion that a 
one-step approach solving a nonsymmetric matrix would be an improvement. 
This would remove the need for under relaxation and much of the data handling 
would be simpler and faster. 

The question of whether to solve a large, sparse, linear system by direct 
means or through relaxation is still very much a controversy in the CFD com- 
munity. We used the direct method and still believe it is the best method on the 
‘use-once’ experimental programs where computational stiffness is anticipated. 
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We did, however, have difficulty with the large cpu time requirement of the 
direct solver and this, more than anything, limited the maximum Reynolds 
number attained. 

The forward-backward solver was very effective as a means of reducing the 
amount of memory required for direct solution of the matrix with only a modest 
increase in cpu time. However the forward-backward solver does have difficulty 

with round-off error amplification in the backsweep. 

* 

The nonlinear scale function used in this study did not enhance convergence 
and in retrospect appears to unnecessarily complicate the problem. 

The Japanese fan adaptive grid was useful in enhancing convergence and 
accuracy. However the stretcher required considerable programmer effort to set 
up and properly tune. No suitable coordinate stretcher in the f direction was 
found that would match the boundary conditions at both zero and infinity and 
give useful control over point placement in the interesting regions of the flow. 
The use of the multipole expansion was found to provide an effective far field 
boundary condition as long as a Neumann boundary condition was used rather 
than a Dirichlet. 

The multipole solution provided some unexpected dividends. One was the 
generation of a quality factor for comparison of the computed value against the 
analytic value of the dipole coefficient. The other dividend was the interesting 
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relationship of multipole coefficients with respect to each other at various Rey- 
nolds numbers. 

A line of speculation developed from this multipole relationship led to the 
belief that a boundary location = 15 is acceptable only for Re < 29.22. 
Otherwise the multipole series is divergent causing the boundary condition to be 
invalid. 

The computer animation proved to be a very useful tool in both understand- 
ing the flows presented and as a teaching aid. The animation revealed features 
about the critical points and their relation to the time lines that could not be 
fully realized using just the entrainment diagrams. 

2. Physics Of Axisymmetric Jets 

One of the main conclusions regarding the physics of impulsively started jets 
and vortex rings involves the discovery of three distinct slates of motion for the 
various flows studied. Transition Reynolds numbers are given in Table 7-1. 
State 1 is for 0 < Re < Re v Application of this flow to an initially 
straight time line leads to a smoothly growing, roughly parabolic front. This 
flow has a single on-axis stable node. State 2 is for Re j < Re < Re 2 . 
Application of this flow to an initially straight time line leads to a curved line 
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Table 7-1 

Transition Reynolds Number 


Flow Type 

First Transition Rej 

Second Transition Re 2 

Stokes Vortex Ring 

18.1749 

23.4105 

Stokes Ramp Jet 

3.7386 

5.7887 

Stokes Round Jet 

6.7804 

10.0909 

Navier Stokes Round Jet 

5.5 ± 0.4 

7.54 ± 0.16 


with a corner located at an angle larger than the angle of the off-axis critical 
point. This flow has an on-axis saddle point and two off-axis stable nodes. State 
3 is for Re 2 < Re. Application of this flow to an initially straight time line 
leads to a geometry that looks like a mushroom or smoke ring with the volutes 
centered on the off-axis critical point. This flow has an on-axis saddle point and 
two off axis stable foci. 

State space analysis was used to define the structure of the flows under 
study. It was shown that employing continuity, boundary conditions and self 
similarity makes it possible to draw the trajectory in p , q space for axisym- 
motric jots, even in the absence of a solution. From a knowledge of this p, q 
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trajectory one can predict the different kinds of flow topologies which will occur 
as the Reynolds number is varied before actually solving the governing equations. 
The p , q trajectory helped us to recognize that State 2 is a distinct state of 
motion in axisymmetric jets. 

Much of the work revolved around the generalized solution of the Stokes 
equation (see App. A). Three creeping solutions (Stokes solutions) were solved 
for and studied in detail. These solutions were examined at higher Reynolds 
number where the creeping approximation is invalid. In fact, state-space analysis 
shows that the Stokes solution of the round jet will have the same topological 
properties as the nonlinear solution even though the approximation is beyond its 
region of validity. There were, however, many differences between thee linear 
and nonlinear round jet solutions. The most. obvious as seen in Table 7-1 is that 
the transition Reynolds numbers are different. Also, the spreading angle behavior 
is completely different. The spreading angle in the Stokes solution goes from 0 
to jt/2 as Re becomes infinite. In the nonlinear solution the spreading angle 
quickly grows with Reynolds number until the second transition. As the Rey- 
nolds number is increased beyond Re 2 the spreading angle reaches a maximum 
of 34.4 degrees and then begins to decrease. For the nonlinear solution the size 
of the vortex increases slowly while the stem of the mushroom quickly grows in 
length. In essence because there is no convective term, the vortices of the Stokes 
solution do not move axially away from the momentum source. 
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Strain in the fluid flow was both an interesting observation and a source of 
difficulty in the computer animation. It was found that the highest strain 
occurred at the on-axis saddle point (see Fig. B-42, Re = 10.7). The strain is 
observable as separation between each of the computed points (represented as 
crosses). Strain decreases as the time line rolls up towards the stable focus from 
the saddle point. The time line is folded at the stable focus and rolled up into a 
vortex centered about the stable focus forming the mushroom shape. The region 
of the time line between the stable focus and the saddle point has more strain 
than the region between the mushrooms and the stable focus. Strain increases 
as the time line rolls up from the stem into the vortex towards the stable focus. 

In the animation of the vortex ring it was found that when the flow changes 
from State 3 to State 2 the vortex stops rolling up. Therefore it is possible in a 
mixing process for the rolling up to terminate before the interface can be 
entrained into the vortex, i.e., it is possible that no mixing occurs at all in a vor- 
tex ring’s vortex. 

One of the most remarkable observations in this study was that there was 
no local accumulation of vorticity at the stable focus of the round jet. The 
greatest vorticity in the round jet occurs in a 1/r 2 singularity at the momen- 
tum source. In the neighborhood of the stable focus the vorticity is smooth and 
decreasing monotonically with radius. This observation is counter to the widely 
accepted belief that vorticity must be accumulate within a vortex. The Stokes 
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ramp jet like the round jet has its vorticity peaked at the momentum source with 
the vorticity smooth through the stable focus. This sort of vorticity behavior will 
also be the case for the other types of axisymmetric jets for integer m > 1 
(Reynolds number increasing with time). For the Stokes vortex ring (m = -1) 
the situation is different. In this flow the vorticity peaks at 
(f,0) =ss (± 1, jt/2) for all Reynolds number. However, the closest a stable 
focus ever gets to this vorticity peak is at Re — ► oo where the focus is at 
(£., d e ) =ss (± 3.0224, jt/2). In the Stokes solution the vorticity peak and criti- 
cal point never occupy the same location. However, it is an open question 
whether the Navier-Stokes vortex ring has a coinciding stable focus and vorticity 
peak. It is quite possible that the inclusion of the convective term may cause the 
vortex to be a site of vorticity accumulation in the vortex ring at high Reynolds 
number. Resolution of this question will have to await future study. 
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Chapter VIII 


SUGGESTIONS FOR FUTURE WORK 


Of the infinite number of flows under the category of axisymmetric jets, one 
flow, the round jet (m = 0), was studied both analytically in the creeping 
approximation and numerically for the full Navier-Stokes solution. Clearly some 
of the other axisymmetric flows are worthy of numerical investigation. 

Recommended flows for further study are the following. 

1. The vortex ring (m = -1); 

2. The ramp jet (m = 1); 

3. Hill’s spherical vortex (m = 3/2); 

4. The two unnamed flows for m = 1/2 and m = - 1/2. 

The spherical vortex should be subjected to a closer analytic investigation 
since it is the author’s opinion that there may exist a whole family of exact 
axisymmetric jet solutions where the convective term goes to zero in the Navier- 
Stokes equation. 

The family of orthonormal solutions (m negative half integer) shown in 
PRECEDING PAGE BLANK NQl FILMED 
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App. A, may also serve as a basis for a spectral investigation of axisymmetric 
jets. 

The nonaxisymmetric round jet is a flow that is very exciting. Like the 
axisymmetric problem the nonaxisymmetric round jet is autonomous and self 
similar. Since this flow is not constrained in the third dimension it can theoreti- 
cally display a vastly more complex set of geometries. This flow would also serve 
as an ideal basis for utilizing three dimensional state space. Three dimensional 
state space is a mathematical tool that will have a much more immediate appli- 
cation to real world flows than the axisymmetric methods developed in this work. 
It will be extremely interesting to see whether the three dimensional state space 
can be used to deduce a trajectory similar to the one discovered in this work 
based only on continuity, boundary conditions, and self similarity. 

The data generated in this study could be further utilized by assuming that 
the fluid is divided into different types with a reaction occurring at the interface. 
The data of this study could drive calculations on the mixing and reaction rates 
of the different fluids. 

There also exists the possibility of analytic solution of the Navier Stokes 
round jet. The inquiries of App. A, coupled with some of the observations of the 
numeric study, represent a set of signposts which indicate that a closed-form 
solution may exist which only awaits the advent of some new mathematical tech- 
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nique to bring about its discovery. 

The problem of computer animation is another area that is very interesting. 
In this study a two-dimensional representation of 1000 points was animated. In 
a three-dimensional study the level of complexity would probably go up by 
several orders of magnitude. Flow visualization actually represents a barrier to 
the advance of understanding in fluid mechanics. For example, the author can 
visualize a two-dimensional flow that is translating and deforming but cannot 
visualize this for a three-dimensional flow. The situation is analogous to a chess 
master who can visualize an entire chessboard and examine scenarios five or six 
moves ahead while a novice player can only visualize three or four pieces and see 
only two moves ahead. For those of us who are not master chess players there 
will be a need for movies showing three-dimensional fluid flows translating and 
deforming so that we can perceive the underlying physics. However, the com- 
puter and peripheral capability for such computer animations does not yet exist 
(the storage, I/O requirements, software are not sufficiently advanced). The jet 
flows described in this work can provide an ideal medium for advancing the tech- 


nology of flow visualization. 
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Appendix A 

JET FUNCTIONS AND THEIR PROPERTIES 


1. Introduction 

Fluid is disturbed by a point momentum source whose time behavior is well 
defined. The flow is constrained to be axisymmetric, and the nonlinear term of the 
Navier-Stokes equation is removed. The equation of motion is the axisymmetric Stokes 
equation which, when expressed in spherical coordinates, is 


d uj 
d t 


v 




(A-l) 


The problem is expressed in spherical polar coordinates (r, 0, <j>) 
where 

9 = 0 is the axis of symmetry, 

t = time 

v = the kinematic viscosity. 
u(r,9,t) = vorticity in the <j> direction, 

Equation (A-l) may be reduced by the introduction of a new dependent variable and a 
self similar coordinate: 


S = 


r 

v/2 ~vt 


(A-2) 


where 


u = Pv- 2 t m - l W(S,0) 


(A-3) 


preceding page blank not filmed 
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S = similarity variable, 

P = characteristic strength of the momentum source and is the 
parameter which, along with u governs the scaling 
properties of the flow 
W = dimensionless vorticity 

A 

m = index derived from the dimensions of P . 

The similarity variable S is a modified version of the earlier similarity variable 
£ used in the main text of this thesis. The two variables are related as follows: 


$ - 5v/2 

Equations (A-2) and (A-3) can be combined with (A-l) to give 


(A-3a) 


d 


+ + 

1 

d 

• a dW 
sm0 "-T- 

d S 

as 

dS 

sin 0 

d 9 

5 0 


(A-4) 


-^-+ 2(m - 1)S 2 
sin 2 0 


W = 0 


Equation (A-4) can be reduced into two ordinary differential equations by separation of 
variables 


j '1 £ i IgA. — \ - 2 

dS* + (5+ dS 2 


2S 2 


J = 0 


(A-4a) 


1 d_ 

sin 0 dO 


sin 0 


dP 

dO 


*('+!)- 


sin 2 0 


P = 0 


(A-4b) 


where 
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W(S,0) = cF(0)JlS) (A-4c) 

and where / and m are real constants, (A-4b) is the Legendre equation of first 
order, a tabulated function. The constant c is for boundary condition matching at 
the far field. For the purposes of this paper we will examine only the solution of (A-4b) 
that corresponds to a dipole flow in the far field, viz., 1 = 1, P(0) = sin0 . we may 
now reexpress (A-4c) as 

W(S,0) = c sin 0 J[S) . (A-5) 

It should be emphasized that J(S) is not a Bessel function (in fact Bessel functions 
never appeared in this study). With 1=1, (A-4a) becomes 


£J 


+ 




■4-+W- 1 

S 2 


7=0. 


(A-6) 


2. DERIVATION OF THE “EASY” SOLUTIONS 

The ordinary differential equation, (A-6) will be the center of study. Solutions of 
this equation will be called “jet functions”. Equation (A-6) by the following substitu- 
tions can be cast into two different forms which are polynomials without transcendental 
function. 


Let 



(A-7) 


giving 
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Let 


giving 


d*E , „ dE 
dS 2 dS 


+ 


1-2 m--2- 

S 2 


E 


0 . 


s* 

US) = ^j~F m (5) 


(A-8) 


(A-9) 


£F s dF 
dS 2 dS 


- 2 



0 . 


(A- 10) 


Equations (A-8) and (A-10) can each be solved by means of the Frobenius method. For 
(A-8) we find 

E{S) = § CjSi +n (A-ll) 

;=o 

where 

n = -1,2 and C x = 0, C 0 7^ 0 ; (A-12a) 


Cj — Cj_ 2 


-j-n + 2m+ 1 

.0 




(A-12b) 


For this O.D.E. there is one solution with a second order zero at the origin, and a 
second solution which is a first order pole. The recursive formula, (A- 12b), can be used 
to find closed form solutions by recognizing that the index j must be an even positive 
integer and that the numerator in (A-12b) goes to zero at some value of j truncating 
the infinite series for a polynomial solution. The results below point to polynomial solu- 


tions: 
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n = -1. first order pole: 


r — r 2(m + 1) — 7 

“ C '- 2 iU- 3) 


(A- 13) 


For polynomials: m = j/2 - 1, therefore m = 0, 1, 2, 3, 4, 


n = 2, second order zero: 


r — r 2m -1- L 

Ci ~ C '- 2 iU+ 3) 


(A- 14) 


For polynomials: m = — [l + j], therefore m 

it 


2 ’ 2 ’2 ’ 2 ’ 


j = 2, 4, 6, 8, 10, ... 


we can go through similar reasoning for (A-10) with the solution being 


F(S) 


E 


j “0 


Cj- s i+n 


where 


n = -1,2 and C x = 0, C 0 7^ 0 , 

r _ r . j± n± 2(m - 1) 

; ; ” 2 y 2 + 2ny-y-2-n+ n 2 


(A- 15) 


(A- 16 a) 
(A-16b) 


This result as before has a solution with a second order zero and a first order pole. The 
polynomial solutions can be recognized in the same manner as before: 


n = -1, first order pole: 
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Ci C ‘- 2 Hi- 3 ) 


(A-17) 


113 5 

For polynomials: m = (3-j)/2, m = — , , 


n = 2, second order zero: 


n — r 2m + j 

C ‘ = C '- 2 Hi + 3 ) 


(A- 18) 


For polynomials: m = - j/2, m = -1, -2, -3, -4, 


j = 2, 4, 6, 8, 10, ... 


These polynomial solutions which shall be called “Easy Solutions” are all expressed 
in closed form as finite series. When the E m and F m polynomials are cast back 
into their original J m form via (A-7) and (A-9), the first order pole solutions of E m 
or F m polynomials become second order poles, which shall be called “B m func- 
tions” and the second order zeros of E m or F m polynomials become first order zeros 
which shall be called “G m solutions”. These independent functions are thus 

= C, G m (5)+ C 2 B m (S) (A- 19) 

where C x and C 2 are constants of integration. 
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2*'- 1 
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1 

(m\)S 


(A-20) 


B 0 (S) = -- 


_1_ 

S 2 


(A-21) 
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where m = 1, 2, 3, 4, ... . 

Rodrigues’ formulae and generating functions have been derived for the 
equations. For (A- 20) 


BAS) = 


- 5 * 
m ~l C 2 


j2m-2 


(2m-2)!mS 2 dS 2m ~ 2 
with generating function ^ where 


(5 2 -2m + 1) e 
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cosh(S\/2X)- 5v/2Xsinh(5v/2X) 


1 

J 


where 


00 

njf)(5,x) = E Bj{sw . 
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(A-22) 

(A-23) 

(A-24) 

(A-25) 

above 

(A-26) 

(A-27) 

(A-28) 


For (A-22) 
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(2m+l)!S 2 dS 2m ~ x 


(S' 2 -2m) e 2 


(A-35) 


n i + J 5 - x > - ik 
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cosh(5\/2\)- — sinh (Sv/2X) (A-36) 

5v2X 


where 


nT’ (5,x) 

T +m 


= £ G Zt AS)\i 


/=o 2 +J 


(A-37) 


where m = 1, 2, 3, 4, ... . 

For the functions described by (A-22) and (A-24) there has been found the follow- 
ing relationships with known special function: 


G^{S) = 


LDll 2 

(2m-l)!\/2 5 2 


4m+2 Him+l 
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where m = 1, 2, 3, 4, ... and H m ^S'/\/2| is a Hermite polynomial and 


I 
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M{-m, 




is a confluent Hypergeometric function. This connection between 


these special functions and jet functions greatly aided the uncovering of the jet 
function’s properties. This particular line of inquiry still promises more interesting 
discoveries for some future work. 


Recursion relations have also been found for the jet function. For 


m < — 
— 9 


it 


was found that 


[2m -5] [m-1] J m _ 2 - [5-4m-S 2 ] J m _ x + 2J m = 0 


(A-42) 


(m < integer or integer) 


For m > 2: 


[2 m -3 ] m J m + [5 -4m -S 2 ] J m _ x 4- 2 J m _ 2 = 0 (A-43) 

(m > 2, integer or — integer) 

2 

J m may be replaced by G m or B m individually in (A-42) and (A-43). 

Other properties of jet functions are more easily examined if we reduce the 
independent functions of G m and B m down to their primitive polynomials. We do 
this by employing (A-7) and (A-9) and defining some new notation. Let 

F„(S) = C,Fi a HS)+ (S) (A-44) 

E m (S) = C t Ei°HS) + C 2 Ei B HS) (A-45) 

where 
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-S> 

a m = Fi c| = -J (A-46) 

-S» 

B m = ■*-£- FW = 1 £<») (A-47) 

and C 1; C 2 are constants of integration as in (A-19). 

Two recursion relations have been found which have a differential component for 

m < -~ 

- 2 

dF m 

[2m -3] m F m _ x = [-2m + 1 - S 2 ] F m + S (A-48) 

Equation (A-48) is true for both F„ and F™ but only if m < -£• . For the 

a 

case of m > 2 the following is true: 

dE m 

-2 £ m _, = [-2m + 1] + 5 — • (A- 49) 

Equation (A-49) is true for both E™ and E ^ but only if m > 2. Imaginary 

relationships analogous to those of the trigonometric functions cos(/0) = cosh(0), 

-»sin(i'0) = sinh(0) have also been found for jet functions (not an unreasonable 

event in light of the generating functions): 


Os) = 

-iF [ B) ( iS ) 

r m 

(A-50) 

^3 

3 

II 

-«• Os) 

(A-51) 
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E ia \ 3 (S) 

7 


f-2-i m 


(A-52) 


F-Zi (S) = E ia) 3 (iS) . 

m + - 


(A-53) 


where m = 0,1, 2, 3, ... 


It was the discovery of these imaginary relationships that largely unlocked the 
study of this problem. Now one need only study the functions of 

m < — (F m polynomials which, because of the alternating sign in the terms mak- 

it 


ing up these functions, were more easily matched to known special functions, (i.e., the 
Hermite polynomial). E m polynomials can be found by simply mapping over from the 
F m polynomials using the above relations. 


3. DERIVATION OF “HARD” SOLUTIONS 

In the previous section we derived what are called “easy solutions”; however, this is 
only half of the story. The reader may have noticed that for m = 0 we have a B 0 
solution but no G 0 solution. Likewise for the rest of the easy solutions, another solu- 
tion needs to be found. These other “hard solutions” as their name implies, are of a 
more complex structure. Hard solutions can be expressed in a single infinite series form 
through application of the Frobenius method. However such a form is virtually useless 
for satisfying boundary conditions or in examining the solution’s properties. The more 
appropriate forms for examining hard solutions are shown below: 
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E m {S) 


DJS) 



-s 1 

e 2 -E m {S)erf 



(A-54) 


F m ( 5 ) 



ii S/v’2 

D m (S) e 2 + F m (5) / e <‘ dx 

0 


(A-55) 


where E m and F m are easy polynomials, E m , F m are hard solutions, and D m 
is a function which will be called a “hard driver polynomial” which will be defined 
shortly. Hard solutions can be converted back to the jet function form J m by using 
(A-7) and (A-9). Easy polynomials of first order zero type (G superscript) always have 
second order pole type (B superscript) hard solutions. Easy polynomials of second 
order pole type (B superscript) always have first order zero type (G superscript) hard 
solutions. B superscript hard solutions have B superscript hard driver polynomials; and 
G superscript hard solutions have G superscript hard driver polynomials. Hard solu- 
tions and hard driver polynomials can be raised to higher or lower order by using the 
recursion equations (A-41) and (A-43) (however no differential recursion formula such 
as (A-49) has been discovered for hard driver functions). The imaginary relations for 
hard solutions are almost the same as for easy solutions: 


F™-1 ( 5 ) = « E m 3 (iS) 

m + - 


e [B) 3 (S) = .• FZ-l (*s) 

m+ 2 




(iS) 


(A-56) 


(A-57) 


(A-58) 
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(S) 



('•?) 


(A-59) 


where m = 0, 1, 2, 3, 4, ... 

The hard driver polynomial D m (S) for (A-54) may be solved by finding the par- 
ticular solution to the following O.D.E., which was derived from inserting (A-54) into 
(A-8): 


££n 

dS 2 


- S 


dD m 

m 

dS 


- 2 


m + 




dE m 

D m = 2 — - 
m dS 


(A-60) 

-s t 


The homogeneous solution to (A-60) is simply F m which when multiplied by e 2 
in (A-54) is E m {S) (the easy solution) and therefore redundant (it can be absorbed 
in the constant of integration). A simple closed form solution for D m (S) is not yet 
attainable. The form presently known is shown below for one case derived from (A-60), 
(m > 0, integer only): 




E C.S "- 2 

Jfc-1 


Cl = 1 

rj _ o m-1 

m (2m-2)! 


C k 


1 

k+ m - 1 


(2k 2 -k-l)C k+l + 


2*(2Ar— l) 2 (m !) 
(2Ar)! (m-k)\ 


(A-61) 


(A-62) 
(A- 63) 


(A- 64) 


where m = 1, 2, 3, ... and C m is used to end the series. The preferred method for 
calculating hard drivers is not to use (A-61) or its counter part for but rather to 
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use the recursion relation (A-43) and use (A-60) to check the result. Below is 
Table A-l of known hard driver functions. 


One may calculate negative m, and m 


1 

2 


hard drivers (for F m ), 


by first calculating the positive m value hard solutions by using the above tech- 
niques and converting them over to their imaginary counterparts via (A-56) and (A- 
59). The hard drivers can be extracted by using the definition of (A-55). However 
for physical reasons, which will be examined in the next section, negative m 
valued hard solutions are of little interest. 


4. INTERLEAVING OF SOLUTIONS 

Table A-2 can now be generated which shows some of the jet function primi- 
tive polynomials and how these functions interleave. The reader should note how 
the easy and hard solutions alternate in both the G m and B m solutions. Also 
note how the J,/ 2 solution (which has a F m polynomial form) penetrates into the 
E m domain, leapfrogging J 0 • An example of how a complete solution is 

extracted from Table A-2 is shown with m = ~ : 

It 

<^5 (S') = C l G^(S)+ C 2 B^(S) 

2 2 2 


(A-65a) 
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Table A-l 

Hard Driver Functions 


D 


(G) 


= -1 


d\° ] = -1 


D[ e) = -i 


5 _ — 

5 5 


Z) 


(G) 


i r 


l 


f i 1 + S2 J 


D? = 
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1 

15 


45 + 5 3 - 




J. 
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1 + f 52 + T s * 
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Interleaving Of Jet Function Solution 
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The imaginary complement of J 5 is J_ 2 which is shown below: 

7 



+ 


~T + 15 1 



dx } . 


(A- 66) | 


At this point boundary conditions for S — ► oo should be examined. It has been 
shown that vorticity from a point disturbance will decay as a Gaussian multiplied by a j 

1 

slowly varying function. In (A-65) we observe that this boundary condition is attain- 
able only if C x = C 2 , which causes the easy solution to disappear at infinity. In 
(A-65b) the solution will converge only if C 2 — 0 . These two cases are representa- 
tive for all J m with polynomial solutions (for the E m solutions C x = C 2 and for the 
F m solutions C 2 = 0 ). It is for this reason that the F m hard solutions are of little 
physical interest in the unbounded problem, for with the exception of the solution 

m = — , , -1 all of F m hard solutions diverge at infinity, and in the cases of 

2 2 

1 1 

2 ’ 2 ’ 


m = 


1 the solutions do not decay as a Gaussian multiplied by a slowly 
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varying function and are therefore inadmissible. Based on this boundary condition we 
can write a more restrictive form for «/ m (S) where we set C l = l for conveni- 
ence: 


J m (S) = — ^ m (5)+ £ m (S)j for m > 0, integer and 1/2 integer ; (A-67) 

-S» 

J m (S) = — for m = ™ < 0, integer and 1/2 integer ; (A- 68) 

S L J - 


5. ORTHOGONALITY PROPERTIES OF JET FUNCTIONS 
It was found that (A-4b) can be recast into the Sturm-Liouville form: 



' £ j,' 


Si Si 

d 

dS 

2 -f J 

+ 

2(1- m)5 2 e 2 - /(/ + 1) e 2 


J 


0 . 


(A- 69) 


This result opens up many possible lines of inquiry, but the only one examined so far is 
that of orthogonality. Orthogonality relations can probably be deduced for all the jet 
functions ( J m ) if one were interested in integrating from zero or infinity to a zero of 
the jet function which had the appropriate boundary condition satisfying the Sturm- 
Liouville theorem. However the only jet functions that are orthogonal over the range 
(0, oo) are the negative m integer type (Fl„). The orthogonality properties of 
these particular functions are defined as follows: 


Let: 
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2 2m (m 2 + m) /T 

(2m +2)! V 2 


From (A-68) 


(A-70) 


= - 



f> f-2V +1 (i 2 + /IS 2 ' 
y=i (m-;)!(2; + 2)! 


for m = 1,2, 3, 4, • • • (A-71) 


We use the weighting function from (A-69). The orthogonality integral is: 


5 s 


/ J_ n (5)S 2 e 2 dS 



m n 
m = n 


(A-72) 


With (A-70), (A-71) and (A-72) as a basis, one could employ the following Fourier form: 


A m = 


/ $(S) J- m {S)S*e 2 dS 


(A-73) 


*(*) = E X 1 - J-mix) (A-74) 

where $(x) is the function to be analyzed. It has not yet been established whether 
(A-74) will pass the completeness test. Assuming that (A-74) is complete, if one had a 
vorticity distribution that can be represented with the self similar coordinate S, one 
could use (A-74) to spectrum analyze the function. This assumes that the resultant 
equations in (A-73) and (A-74) are convergent, which will probably be the case only if 
the Reynolds number of the function under consideration ^(x) is decaying in time 
since (as we shall see in the next section) all of the basis functions J_ m (S) employed 
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in (A-73) are of this Reynolds number behavior. 

6. STREAM FUNCTIONS, VELOCITIES AND REYNOLDS NUMBERS 

FROM JET FUNCTIONS 

With the method of generating the vorticity functions established, the means 
by which to calculate the stream function, velocity, and Reynolds number, along 
with deducing the forcing function should also be examined. Following the method 
used by Cantwell [1980] the velocities and stream function are converted to dimen- 
sionless quantities: 

u = Pv-W t m ~ 1 / 2 U{S,9) (A-75) 

v = Pv~ 3 / 2 t m ~ 1 / 2 V{S,6) (A-76) 

ip = Pt/-i/2 r +1 / 2 g{S,6) (A-77) 

where u is the radial velocity with U its dimensionless counterpart, v the 
velocity in the theta direction with V its dimensionless counterpart and the rest 
of the variables the same as used in (A-3). The definition for stream function and 
vorticity in axisymmetric spherical coordinates is: 
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1 

8 tp 

(A-78) 

i^sin 6 

8 9 

-1 

8 rp 

(A-79) 

r sin 0 

8 r 

-k < re > 

8 u 

~ d 9 ' 

(A-80) 


The previous six equations are used with (A-2) and (A-3) giving 


-2 3 / 2 S W = 


1 d 2 g | 1 d 

sin 9 d 5 s & 89 


1 d g 
sin 9 8 9 


(A-81) 


It can be shown that an axisymmetric stream function for a Stokes vorticity solution 
that is consistent with (A-5) (/ = 1, dipole case), must be of the following form: 


g(S,9) = c sin 2 6 R(S) . (A-82) 

Combining (A-83) with (A-81) and (A-5) we find that 


2 3/j AS) = ~ 


i Ts < 5 *> 


(A-83) 


where J(S) is the jet function which we just solved. Equation (A-83) is readily 
integrated, yielding 


fl(S| 


= + Ci & - [ /<?( J 4S)dS) rfs] 


(A-84) 


where C x and C 2 are constants of integration, governing the irrotational dipole com 
ponents of the stream function. The dimensionless stream function, and velocities can 
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now be expressed utilizing (A-75), (A-76), (A-82) and (A-84): 


g(S,9) = sin 2 0 t —7 + C 2 5 2 - 




— [ /5 2 ( / ./(S)^)^] | .(A-85) 


U = 


V = 


1 

<L± 

25 2 sin 6 

d 9 

-1 

<L± 

25 sin 9 

d S 


(A-86) 


(A-87) 


We may derive the Reynolds number from the unsteady particle paths given by 


dr 

dt 


= u 


(A-88) 


d± 

dt 


v_ 

r 


A logarithmic relation for time is employed: 

7 = In t 

Also, a Reynolds number will be used of the following form: 


(A-89) 


(A-90) 


Re = 


Pt m 


II 

2 


(A-91) 


Equations (A-75), (A-76) and (A-2) are used with (A-88), (A-89), (A-90) and (A-91) giv- 
ing 


— = Re 2 U - — 

dr 2 


(A-92) 
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d0_ 

dr 



(A-93) 


The forcing function at the point momentum source of the jet can be deduced by 
first recognizing that the impulse of the jet is conserved. Since impulse is conserved we 
may write down the impulse integral as an equation of constraint: 


where 


I_ = 

P 

is impulse, and 


7TOO 

I II <« cos 0 - v sin 0) 2m 2 sin 0 dr d0 (A-94) 

2 o o 

p is density. The definition of impulse is the following: 


I_ 

P 


- J Fdt 
P o 


where F is the force applied to the fluid at the momentum source. 
We combine (A-94), (A-95) with (A-75), (A-76) and (A-2) giving 


(A-95) 


l r~ 7T OO 

1 A O./O A 

- f F dt = P r +1 // {U cos 0 - V sin 0) 2it& sin 0 dS dO 

P 0 2 0 0 

(A-96) 


If the components of (A-96) that are not functions of time are collected one can define 
the following constant: 


K 


3\/2 

2 


7TOO 

p P f f ( U cos 0 - V sin 0)2x5? sin# dSd0 

o o 


(A-97) 


Equation (A-96) can now be reexpressed as 
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t 

f jr jf — l i a not 

j * ***> — * • yr\-vo j 

0 

With (A-98) we may deduce the force function. Let us examine some examples. 

Let: 


F = Pt . (A-99) 

We find that (A-98) is satisfied only if m = 1. For the case m = 1 the Reynolds 
number described by (A-91) increases proportionally to f 1 / 2 . 

Let: 


F = Pu{t) (A-100) 

where u[t) is a Heaviside step function. For this case m = 0, the Reynolds 
number is not a function of time since time is to the zeroth power. The case of m = 0 
is unique in that it is autonomous. This solution corresponds to a round jet. Let: 

F = P6{t) (A-101) 

where <!>(<) is the Dirac delta function. For this case m = -1 since time is removed 
from the left-hand side of (A-99). The Reynolds number decreases proportionally to 
l/f 1 / 2 making this case complementary to the case of m = 1, (A-99). Let: 

F = P l 3 / 2 (A- 102) 

3 . 

For this case m = 3/2 , the Reynolds number increases proportionally to t 4 . 
This solution corresponds to the “Hill’s spherical vortex”. Let: 
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F = Pt 1 / 2 (A- 103) 

j. 

For this case m = 1/2 , the Reynolds number increases proportionally to t 4 . 
This solution is the imaginary complement of the m = 0, round jet solution. Unlike 
the round jet, the “m = 1/2 jet” has never been extensively investigated. Let: 

F = P r 1 / 2 (A- 104) 

TTr«r fine o ooo m — _1 /9 t.hp R pvrmlde: numhpr deerpases nmnortinnallv t.n 1 It. 4 / 4 

This particular solution is very unusual with regards to its Stokes solution stream func- 
tion which appears to be the simplest of all the dipole solutions examined. This solution 
has never been extensively investigated. 

It should be noted that if one integrates in time (A-101) (the case for m = -1), 
then (A-100) (case m = 0), is obtained. Likewise if (A-100) is integrated in time, 
then (A-99), (case m = 1) appears. We also find that (A-104) integrated yields (A- 
103) and (A-103) integrated gives (A-102). It should also be noted that with the excep- 
tion of m = 0, for negative m the Reynolds number decays with time, while for 
positive m the Reynolds number increases with time. The round jet (m = 0) does 
neither because it is autonomous. The case of m = 1/2 has an increasing Reynolds 
number but it is of an F m polynomial type (see Table A-2). This type is of decaying 
Reynolds number with this single exception. Needless to say, m = 1/2 is a very 


strange solution indeed. 
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7. JET FUNCTIONS AND THE NAVIER STOKES EQUATIONS 

In the previous section it was observed that in the six examples shown for 
m = 3/2, 1, 1/2, 0, -1/2, -1 three are of classical interest; i.e., the round jet, 

vortex ring, and spherical vortex. It should be recalled that all of the solutions exam- 
ined so far have been for the dipole case, / = 1. This document would be too long if 
the cases of 1=2, 3, 4, ... were also examined thoroughly, but for the sake of 
interest let us examine the family of easy functions of the form: 

J = S Q (A- 105) 

where o is an integer. Equation (A-105) is inserted into (A-4a) yielding 

a(a + 1) - /(/+ 1) = 0 (A- 106) 

a + 2 = 2m . (A-107) 

From (A-106) and (A-107) we derive Table A-3. 

It is interesting to note that the two solutions admitted for / = 1 in Table A-3 
both have analytic solutions to the Navier Stokes equations. Hills’ spherical vortex 
satisfies not only the Navier Stokes equation but the Stokes equation as well because the 
convective term goes to zero identically for all Reynolds number. The method by which 
the spherical vortex solution is found starts with the following equation for the convec- 
tive term: 


d f , L d r i 

T7l r “ a '] + 77 M = 


(A- 108) 


Equations (A-78) and (A-79) are included with (A-108) giving 
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d 

u> 8 rp 

n r „ . i 

a u> o if) 

d r 

r sin 0 d 0 

6 0 r sin 0 dr 


(A- 109) 


Any Stokes solution that satisfies (A-109) will also satisfy the Navier Stokes equations. 
The vorticity solution for M = 3/2 has been derived in the previous section (Table 
A-2) and is 


Inserting 


w^(r, 0 ) 

2 

(A- 110) into (A-109) gives 


cP 

\/l8 i/ 5 


r sin 0 . 


the relation for the stream function 


(A- 110) 


6 ip 

6 

6 tp 

TJ 

6 0 

6 r 


0 (valid only for m = 3/2) . (A-lll) 


From (A-lll) we find that any stream function consistent with (A-110), using (A-84) 
and the entire family of irrotational axisymmetric stream functions, will satisfy the 
Navier Stokes equation for this case. This complete solution is 


ip(r,0) 


V * potential "b 


-2 1 / 2 c P r 4 sin 2 6 
60 ifil 2 


(A- 112) 


where ip potential is the family of irrotational axisymmetric stream functions: 


^ potential t) = E [ H } (t) S-i + Q j (t) S>' +i ] sin 0 P) (cos 0) (A-113) 

/- 1 L J 

where H } and Q } - are time dependent coefficients determined by boundary condi- 
tions, and Pj is the first order Legendre polynomial. 
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While the solutions of (A-110) and (A-112) have been known for almost a century 
the fact that higher order irrotational stream functions (quadrupole, hexapole, etc.), 
were also solutions seems to have been overlooked. The classical references (Lamb, 
Milne, etc) make no mention of this and Turner [1964] who attempted experimental 
verification of the spherical vortex used only the dipole irrotational component even 
though his experimental data was indicating higher order terms by producing an oblate 
geometry. Also no one appears to be aware of the fact that the spherical vortex can be 
produced by a point momentum source of strength proportional to t 6 !' £ and a Rey- 
nolds number proportional to f 3 / 4 . 

Hills’ spherical vortex is a member of a larger body of solutions (where (A- 109) is 
true), which satisfies both the axisymmetric Stokes equation and Navier stokes equation 
for all Reynolds number, thus being self similar with respect to S, (A-2) for all Rey- 
nolds number. All axisymmetric jet solutions will be self similar with S in the creep- 
ing limit but only one other family of solutions will have this self similarity for all Rey- 
nolds number. This other family is made up of the autonomous solutions like the round 
jet, radial jet, etc., which have been described by Cantwell [1981]. The convective 
term in the autonomous solutions does not go to zero. Even though these solutions are 
self similar with S they satisfy the Stokes equation only in the creeping limit. The 
complement to the spherical vortex in Table 3, the round jet is the creeping limit of 
the Navier Stokes solution discovered by Landau and independently by Squire and as 
mentioned before is an autonomous solution. The Landau Squire solution is generally 
represented as follows: 
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u 


vr 


2 sin 2 0 
A - cos 9 


4WA 2 -l)sin 9 
r 2 (A - cos 9 ) 3 


(A- 11 4) 


(A-115) 


I 


Re 2 
16 x 


A + 


1 A 
3 A 2 -l 



1a+i ' 
U ’ 1 J 


(A- 11 6) 


In working with the first order Legendre polynomial it was observed that the gen- 
erating function for P l l {co$9) , (which was found in Morse and Feshbach [1953]) is 
similar to (A-115). The generating function is 




I 

! 

i 


I 


sin 9 

1 

[1+ h 2 -2h cos0] 2 

Which ever series converges for the given h is the appropriate form. Equation (A- 
117) is very intriguing not only because it represents a special case of the irrotational 
solution (A-113 with constant coefficients), but also because it so closely resembles (A- 
115). This aspect, coupled with the fact that the complete solution of the spherical vor- 
tex stream function included all of first order Legendre polynomials, lead to the desire 
to recast (A-114) as a series of Legendre polynomials. After considerable algebra this 
was done giving the following: 


= = E Pi * i («» *) 


n =0 


n *=0 ^ 

k j small | h | large 


(A- 117) 
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*(r,t) = ^ E ^ [a(4/+3)C„- +1 PJ /+1 <«*»> 


+ (4; -h 5) C 2 y + 2 ^2; + 2 ( cos 


(A- 118) 


c = ^ (2fn + j 1)! (n+ 2 + 2j)! {jn + 1) 

2j+1 „_o (»+ })'■ (2[n+ 2+ 2;])! A 2 " 


(A-119) 


“ (2 [n +J_± 11)! (n + 3+ 2? )! {jn + 1) 

n=0 (n + j + 1)! (2[n + 3 + 2>])! A 2 * 


(A- 120) 


In many ways (A-118) is a more interesting form than its generating function form in 
(A-114) because the solution’s components can be picked out individually and examined. 
It is also interesting to note that the solution of (A-118) has split into odd poles and 
even poles. This leads to a line of inquiry as to whether any of the components will 
independently satisfy the Navier Stokes equations. This question still remains 
unanswered. The lowest order solutions of (A-119) and (A-120) was cast into an alge- 
braic form: 




(A-121) 


A* J_ 

T Cl '3A 


(A-122) 


Equation (A-121) corresponds with the dipole solution and goes to — — as Reynolds 
number goes to zero. Equation (A-122) corresponds with the quadrupole solution and 
disappears as C j approaches — — . It is interesting to note the similarity between 
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(A-121) and (A-116) whether any new insights can be gained from this similarity 
remains to be determined. 

8. CONCLUSIONS AND FUTURE LINES OF INQUIRY 

The methods described in this document will generate low Reynolds number solu- 
tions for any axisymmetric jet. These methods also provide a means of classifying 
different flow types, showing their interrelationships, and isolating particularly interest- 
ing flows. Many of these flows have never been investigated and are each worthy of 
investigation of their unique properties. Perhaps the most interesting aspect of the jet 
functions is their utility in devising solutions to the Navier Stokes equations. There 
exists a family of solutions of which one member is known (Hills’ spherical vortex), 
where Stokes solutions will also satisfy the Navier Stokes equations. Table 3 might 
serve as a sign post in discovering the other members of this family. There also exists 
the family of autonomous solutions of which the round jet is a member. It is possible 
that other solutions can be extracted from the round jet solution cast as an infinite 
series of Legendre polynomials. Higher order components and Reynolds number 
behavior was observed in the spherical vortex. Both the spherical vortex and the round 
jet were observed to be composed of an infinite number of Legendre polynomials. Three 
major payoffs that may come from further investigation in jet functions are: 

1. A generalized solution of the family of flows where the convective term goes 
to zero identically as is the case of the spherical vortex. 

2. Discovery of a far field (hard function) solution to the Navier Stokes equa- 
tion. All solutions presently known are near field, easy functions. 
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3. Development of spectral methods from first order Legendre polynomials 
and/or jet functions for use in numeric computation of unbounded flows. 



Appendix B 


PLOTS AND FIGURES 


This chapter provides the graphical presentation of the results of this study. 



VORTICITY CONTOUR PLOT FOR Re= 0.1 

30 X 30 MESH FOR f.=150 



Fig. B-la Re = 0.1 For 30 X 30 Mesh 
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Fig. B-lb. Re = 1.0 For 30 X 30 Mesh at = 15.0. 











STREAM FUNCTION CONTOUR PLOT FOR Re= 3.0 I VORTICITY CONTOUR PLOT FOR Re- 3.0 

30 X 30 MESH FOR {.=15.0 30 x 30 MESH F0R L =15 ° 


ORIGINAL PAGE fS 

QF_ -PQQR_QUAL1T^ 



Fig. B-3. Re = 3.0 For 30 X 30 Mesh at = 15.0. 






STREAM FUNCTION CONTOUR PLOT FOR Re= 4.0 VORTICITY CONTOUR PLOT FOR Re 

30 X 30 MESH FOR (.=15.0 30 X 30 MESH FOR (.=15.0 







STREAM FUNCTION CONTOUR PLOT FOR Re= 5.0 VORTICITY CONTOUR PLOT FOR Re= 5.0 

30 X 30 MESH FOR £_=15.0 30 X 30 MESH FOR £„=15.0 



Fig. B-5. Re = 5.0 For 30 X 30 Mesh at 
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STREAM FUNCTION CONTOUR PLOT FOR Re= 8.0 VORTICITY CONTOUR PLOT FOR Re= 8.0 

30 X 30 MESH FOR f_=l5.0 30 X 30 MESH FOR f_=15.0 



Fig. B-8. Re = 8.0 For 30 X 30 Mesh at ^ = 15.0 
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STREAM FUNCTION CONTOUR PLOT FOR Re=10.0 I VORTICITY CONTOUR PLOT FOR Re-10.0 

30 X 30 MESH FOR £.=15.0 30 x 30 MESH F0R ^- =15 ° 



Fig. B-10. Re = 10.0 For 30 X 30 Mesh at = 15.0 




STREAM FUNCTION CONTOUR PLOT FOR Re=11.0 VORTICITY CONTOUR PLOT FOR Re=11.0 

30 X 30 MESH FOR f.=15.0 30 X 30 MESH FOR f.=150 
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Fig. B-ll. Re = 11.0 For 30 X 30 Mesh at = 15.0 














STREAM FUNCTION CONTOUR PLOT FOR Re=15.0 VORTICITY CONTOUR PLOT FOR Re=15.0 VORTICITY*(PSI‘*3) PLOT FOR Re=15.0 

30 X 30 MESH FOR £.=15.0 30 X 30 MESH FOR £.=15.0 PLOTTED FROM 30 X 30 MESH FOR PS1NF=15.0 
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Fig. B-15. Re = 15.0 For 30 X 30 Mesh at ^ = 15.0. 




Fig. B-16. Re = 16.0 For 30 X 30 Mesh at = 15.0. 


STREAM FUNCTION CONTOUR PLOT FOR Re=I7.0 VORTICITY CONTOUR PLOT FOR Re=17.0 V0RTICITY*(PSl**3) PLOT FOR Re=17.0 

30 X 30 MESH FOR ?_=15.0 30 X 30 MESH FOR £_=15.0 PLOTTED FROM 30 X 30 MESH FOR PSINF=15.0 


189 


ORIGINAL PAGE TS 
OF POOR QUALITY 



Fig. B-17. Re = 17.0 For 30 X 30 Mesh at = 15.0. 


STREAM FUNCTION CONTOUR PLOT FOR Re=18.0 VORTICITY CONTOUR PLOT FOR Re=18.0 V0RTICITY*(PS1 M 3) PLOT FOR Re=18.0 

30 X 30 MESH FOR {.=150 30 X 30 MESH FOR f.=150 PLOTTED FROM 30 X 30 MESH FOR PS1NF=15.0 
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Fig. B-18. Re = 18.0 For 30 X 30 Mesh at ^ = 15.0. 


STREAM FUNCTION CONTOUR PLOT FOR Re=19.0 VORTICITY CONTOUR PLOT FOR Re=19.0 VORTICITY*(PSI**3) PLOT FOR Re=19.0 

30 X 30 MESH FOR £=15.0 30 X 30 MESH FOR $.=15.0 PLOTTED FROM 30 X 30 MESH FOR PSINF=15.0 
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Fig. B-19. Re = 19.0 For 30 X 30 Mesh at = 15.0. 



Fig. B-20. Re = 20.0 For 30 X 30 Mesh at = 15.0 



STREAM FUNCTION CONTOUR PLOT FOR Re=21.0 VORTICITY CONTOUR PLOT FOR Re=21.0 VORTICITY*(PSI*‘3) PLOT FOR Re=2L0 

30 X 30 MESH FOR {.=15.0 30 X 30 MESH FOR (,=15.0 PLOTTED FROM 30 X 30 MESH FOR PSINF=I5.0 



Fig. B-21. Re = 21.0 For 30 X 30 Mesh at = 15.0. 





STREAM FUNCTION CONTOUR PLOT FOR Re=22.0 VORTICITY CONTOUR PLOT FOR Re=22.0 VORTICITY*(PSl**3) PLOT FOR Re=22.0 

30 X 30 MESH FOR f =15.0 30 X 30 MESH FOR ^=15.0 PLOTTED FROM 30 X 30 MESH FOR PSINF=15.0 
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Fig. B-22. Re = 22.0 Fop 30 X 30 Mesh at = 15.0. 



Fig. B-23. Re = 23.0 For 30 X 30 Mesh at £*» = 15.0. 
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Fig. B-26. Re = 26.0 For 30 X 30 Mesh at £oo = 15.0 





Fig. B-27. Re = 27.0 For 30 X 30 Mesh at 





3) PLOT FOR Re=28.0 



Fig. B-28. Re = 28.0 For 30 X 30 Mesh at 











Fig. B-30. Re = 30.0 For 30 X 30 Mesh at 








STREAM FUNCTION CONTOUR PLOT FOR Re= 4.0 VORTICITY CONTOUR PLOT FOR Re= 4.0 

60 X 60 MESH FOR ?_=15.0 60 X 60 MESH FOR £_=15.0 
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Fig. B-31. Re = 4.0 For 60 X 60 Mesh at = 15.0. 








STREAM FUNCTION CONTOUR PLOT FOR Re= 6.0 VORTICITY CONTOUR PLOT FOR Re= 6.0 

60 X 60 MESH FOR f =15.0 60 X 60 MESH FOR £_=15.0 
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STREAM FUNCTION CONTOUR PLOT FOR Re=10.0 VORTICITY CONTOUR PLOT FOR Re=10.0 

60 X 60 MESH FOR £_=15.0 60 X 60 MESH FOR |_=15.0 
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Fig. B-33. Re = 10.0 For 60 X 60 Mesh at = 15.0. 





Fig. B-34. Re = 15.0 For 60 X 60 Mesh at foo = 15 0. 






ORIGINAL PAGE IS 
OF POOR QUALITY 






Fig. B-35. Re = 17.0 For 60 X 60 Mesh at Coo = 15 0 







STREAM FUNCTION CONTOUR PLOT FOR Re=23.0 VORTICITY CONTOUR PLOT FOR Re=23.0 V0RTICITY‘(PSI*‘3) PLOT FOR Re=23.0 

60 X 60 MESH FOR £.=15.0 60 X 60 MESH FOR £.=15,0 PLOTTED FROM 60 X 60 MESH FOR PS!NF=150 



Fig. B-37. Re = 23.0 For 60 X 60 Mesh at = 15.0. 
















Fig. B-39. Re = 27.0 Fop 60 X 60 Mesh at ^ = 15.0. 




Fig. B-40. Re = 30.0 Fop 60 X 60 Mesh at 
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Fig. B-41. Comparison Of The Stable Node Critical Point With the Saddle 
Point. 
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Appendix C 


SOFTWARE 


Appendix C contains listings of the following software programs. 

• MAVIN 

• INVORT 

• HAMMER 

• PLOT 

• CONTOUR 

• CONVORT 

• VORT3 

• MOVIE 

• RING 

• RING (Script Program). 

These programs are contained in Ref. 13. If Appendix C is desired, please 
contact Professor Brian Cantwell, Dept. Aeronautics and Astronautics, Stanford 
University, Stanford, CA. 94305. 
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